From 23c6f8087e88818e50957afc1a71e20a38c0fc2b Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Thu, 3 Oct 2024 17:07:28 +0000 Subject: [PATCH] build based on 3d167b4 --- dev/.documenter-siteinfo.json | 2 +- dev/api/index.html | 84 +++- dev/index.html | 2 +- dev/objects.inv | Bin 530 -> 535 bytes dev/search_index.js | 2 +- dev/usage/{86e74be1.svg => bb1f32c8.svg} | 520 +++++++++++------------ dev/usage/index.html | 26 +- 7 files changed, 348 insertions(+), 288 deletions(-) rename dev/usage/{86e74be1.svg => bb1f32c8.svg} (54%) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index ac46fbb..cc50bec 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-10-03T15:27:55","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.5","generation_timestamp":"2024-10-03T17:07:25","documenter_version":"1.7.0"}} \ No newline at end of file diff --git a/dev/api/index.html b/dev/api/index.html index c3b092b..b8981e6 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,51 +1,91 @@ -API · SmoQyKPMCore.jl
GitHub

API

Note that for many methods there are two verions, one that relies on taking an instance of the KPMExpansion type as an argument, and a lower level one that does not.

SmoQyKPMCore.KPMExpansion — Type
mutable struct KPMExpansion{T<:AbstractFloat, Tfft<:FFTW.r2rFFTWPlan}

A type to represent the Chebyshev polynomial expansion used in the KPM algorithm.

Fields

  • M::Int: Order of the Chebyshev polynomial expansion.
  • bounds::NTuple{2,T}: Bounds on eigenspectrum in the KPM algorithm.
  • buf::Vector{T}: The first M elements of this vector of the Chebyshev expansion coefficients.
  • r2rplan::Tfft: Plan for performing DCT to efficiently calculate the Chebyshev expansion coefficients via Chebyshev-Gauss quadrature.
source
SmoQyKPMCore.KPMExpansion — Method
KPMExpansion(func::Function, bounds, M::Int, N::Int = 2*M)

Initialize an instance of the KPMExpansion type to approximate the univariate function func, called as func(x), with a order M Chebyshev polynomial expansion on the interval bounds[1] < x bounds[2]. Here, N ≥ M is the number of points at which func is evaluated on that specified interval, which are then used to calculate the expansion coeffiencents via Chebyshev-Gauss quadrature.

source
SmoQyKPMCore.kpm_update! — Function
kpm_update!(
+API · SmoQyKPMCore.jl
GitHub
API
Note that for many methods there are two verions, one that relies on taking an instance of the KPMExpansion type as an argument, and a lower level one that does not.
SmoQyKPMCore.KPMExpansion — Type
mutable struct KPMExpansion{T<:AbstractFloat, Tfft<:FFTW.r2rFFTWPlan}
A type to represent the Chebyshev polynomial expansion used in the KPM algorithm.
Fields
  • M::Int: Order of the Chebyshev polynomial expansion.
  • bounds::NTuple{2,T}: Bounds on eigenspectrum in the KPM algorithm.
  • buf::Vector{T}: The first M elements of this vector of the Chebyshev expansion coefficients.
  • r2rplan::Tfft: Plan for performing DCT to efficiently calculate the Chebyshev expansion coefficients via Chebyshev-Gauss quadrature.
source
SmoQyKPMCore.KPMExpansion — Method
KPMExpansion(func::Function, bounds, M::Int, N::Int = 2*M)
Initialize an instance of the KPMExpansion type to approximate the univariate function func, called as func(x), with a order M Chebyshev polynomial expansion on the interval bounds[1] < x bounds[2]. Here, N ≥ M is the number of points at which func is evaluated on that specified interval, which are then used to calculate the expansion coeffiencents via Chebyshev-Gauss quadrature.
source
SmoQyKPMCore.kpm_update! — Function
kpm_update!(
     kpm_expansion::KPMExpansion{T}, func::Function, bounds, M::Int, N::Int = 2*M
-) where {T<:AbstractFloat}
In-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, expansion order M, and number of points at which the expanded function is evaluated when computing the expansion coefficients. This includes recomputing the expansion coefficients.
source
SmoQyKPMCore.kpm_update_bounds! — Function
kpm_update_bounds!(
+) where {T<:AbstractFloat}
In-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, expansion order M, and number of points at which the expanded function is evaluated when computing the expansion coefficients. This includes recomputing the expansion coefficients.
source
SmoQyKPMCore.kpm_update_bounds! — Function
kpm_update_bounds!(
     kpm_expansion::KPMExpansion{T}, func::Function, bounds
-) where {T<:AbstractFloat}
In-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, recomputing the expansion coefficients.
source
SmoQyKPMCore.kpm_update_order! — Function
kpm_update_order!(
     kpm_expansion::KPMExpansion, func::Function, M::Int, N::Int = 2*M
-)
In-place update the expansion order M for an instance of KPMExpansion, recomputing the expansion coefficients. It is also possible to udpate the number of point N the function func is evaluated at to calculate the expansion coefficients.
source
SmoQyKPMCore.kpm_coefs — Function
kpm_coefs(func::Function, bounds, M::Int, N::Int = 2*M)
Calculate and return the Chebyshev expansion coefficients. Refer to kpm_coefs! for more information.
source
SmoQyKPMCore.kpm_coefs! — Function
kpm_coefs!(
+)
In-place update the expansion order M for an instance of KPMExpansion, recomputing the expansion coefficients. It is also possible to udpate the number of point N the function func is evaluated at to calculate the expansion coefficients.
source
SmoQyKPMCore.kpm_coefs — Function
kpm_coefs(func::Function, bounds, M::Int, N::Int = 2*M)
Calculate and return the Chebyshev expansion coefficients. Refer to kpm_coefs! for more information.
source
SmoQyKPMCore.kpm_coefs! — Function
kpm_coefs!(
     coefs::AbstractVector{T}, func::Function, bounds,
     buf::Vector{T} = zeros(T, 2*length(coefs)),
     r2rplan = FFTW.plan_r2r!(zeros(T, 2*length(coefs)), FFTW.REDFT10)
-) where {T<:AbstractFloat}
Calculate and record the Chebyshev polynomial expansion coefficients to order M in the vector ceofs for the function func on the interval (bounds[1], bounds[2]). Let length(buf) be the number of evenly spaced points on the interval for which func is evaluated when performing Chebyshev-Gauss quadrature to compute the Chebyshev polynomial expansion coefficients.
source
SmoQyKPMCore.kpm_moments — Function
kpm_moments(
-    M::Int, A, bounds, R::T, tmp = zeros(eltype(R), size(R)..., 3)
-) where {T<:AbstractVecOrMat}
Calculate and return the first $M$ moments
\[\mu_m = \langle R | T_m(A) | R \rangle\]
if $R$ is a vector, and if $R$ is a matrix then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A) | R_n \rangle,\]
where $| R_n \rangle$ is the n'th column of $R$.
source
kpm_moments(
-    M::Int, A, bounds, U::T, V::T, tmp = zeros(eltype(V), size(V)..., 3)
-) where {T<:AbstractVecOrMat}
Calculate and return the first $M$ moments
\[\mu_m = \langle U | T_m(A) | V \rangle\]
if $U$ and $V$ are vector, and if they are matrices then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A) | V_n \rangle,\]
where $| U_n \rangle$ and  $| V_n \rangle$ are are the n'th columns of each matrix.
source
SmoQyKPMCore.kpm_moments! — Function
kpm_moments!(
-    μ::AbstractVector, A, bounds, R::T, tmp = zeros(eltype(R), size(R)..., 3)
-) where {T<:AbstractVecOrMat}
Calculate the moments
\[\mu_m = \langle R | T_m(A) | R \rangle\]
if $R$ is a vector, and if $R$ is a matrix then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A) | R_n \rangle,\]
where $| R_n \rangle$ is the n'th column of $R$.
source
kpm_moments!(
-    μ::AbstractVector, A, bounds, U::T, V::T, tmp = zeros(eltype(V), size(V)..., 3)
-) where {T<:AbstractVecOrMat}
Calculate the moments
\[\mu_m = \langle U | T_m(A) | V \rangle\]
if $U$ and $V$ are vector, and if they are matrices then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A) | V_n \rangle,\]
where $| U_n \rangle$ and  $| V_n \rangle$ are are the n'th columns of each matrix.
source
SmoQyKPMCore.kpm_density — Function
kpm_density(
+) where {T<:AbstractFloat}
Calculate and record the Chebyshev polynomial expansion coefficients to order M in the vector ceofs for the function func on the interval (bounds[1], bounds[2]). Let length(buf) be the number of evenly spaced points on the interval for which func is evaluated when performing Chebyshev-Gauss quadrature to compute the Chebyshev polynomial expansion coefficients.
source
SmoQyKPMCore.kpm_moments — Function
kpm_moments(
+    M::Int, A, bounds, R::T,
+    tmp = zeros(eltype(R), size(R)..., 3)
+) where {T<:AbstractVecOrMat}
Calculate and return the first $M$ moments
\[\mu_m = \langle R | T_m(A) | R \rangle\]
if $R$ is a vector, and if $R$ is a matrix then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A) | R_n \rangle,\]
where $| R_n \rangle$ is the n'th column of $R$.
source
kpm_moments(
+    M::Int, A, bounds, U::T, V::T,
+    tmp = zeros(eltype(V), size(V)..., 3)
+) where {T<:AbstractVecOrMat}
Calculate and return the first $M$ moments
\[\mu_m = \langle U | T_m(A) | V \rangle\]
if $U$ and $V$ are vector, and if they are matrices then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A) | V_n \rangle,\]
where $| U_n \rangle$ and  $| V_n \rangle$ are are the n'th columns of each matrix.
source
SmoQyKPMCore.kpm_moments! — Function
kpm_moments!(
+    μ::AbstractVector, A, bounds, R::T,
+    tmp = zeros(eltype(R), size(R)..., 3)
+) where {T<:AbstractVecOrMat}
Calculate the moments
\[\mu_m = \langle R | T_m(A) | R \rangle\]
if $R$ is a vector, and if $R$ is a matrix then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle R_n | T_m(A) | R_n \rangle,\]
where $| R_n \rangle$ is the n'th column of $R$.
source
kpm_moments!(
+    μ::AbstractVector, A, bounds, U::T, V::T,
+    tmp = zeros(eltype(V), size(V)..., 3)
+) where {T<:AbstractVecOrMat}
Calculate the moments
\[\mu_m = \langle U | T_m(A) | V \rangle\]
if $U$ and $V$ are vector, and if they are matrices then
\[\mu_m = \frac{1}{N} \sum_{n=1}^N \langle U_n | T_m(A) | V_n \rangle,\]
where $| U_n \rangle$ and  $| V_n \rangle$ are are the n'th columns of each matrix.
source
SmoQyKPMCore.kpm_density — Function
kpm_density(
     N::Int, μ::T, W
-) where {T<:AbstractVector}
Given the KPM moments $\mu$, evaluate the corresponding spectral density $\rho(\epsilon)$ at $N$ energy points $\epsilon$, where $W$ is the spectral bandwidth or the bounds of the spectrum. This function returns both $\rho(\epsilon)$ and $\epsilon$.
source
SmoQyKPMCore.kpm_density! — Function
kpm_density!(
+) where {T<:AbstractVector}
Given the KPM moments $\mu$, evaluate the corresponding spectral density $\rho(\epsilon)$ at $N$ energy points $\epsilon$, where $W$ is the spectral bandwidth or the bounds of the spectrum. This function returns both $\rho(\epsilon)$ and $\epsilon$.
source
SmoQyKPMCore.kpm_density! — Function
kpm_density!(
     ρ::AbstractVector{T},
     ϵ::AbstractVector{T},
     μ::AbstractVector{T},
     bounds,
     r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01)
-) where {T<:AbstractFloat}
Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here bounds is bounds the eigenspecturm, such that the bandwidth is given by W = bounds[2] - bounds[1].
source
kpm_density!(
+) where {T<:AbstractFloat}
Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here bounds is bounds the eigenspecturm, such that the bandwidth is given by W = bounds[2] - bounds[1].
source
kpm_density!(
     ρ::AbstractVector{T},
     ϵ::AbstractVector{T},
     μ::AbstractVector{T},
     W::T,
     r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01),
-) where {T<:AbstractFloat}
Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here $W$ is the spectral bandwidth.
source
SmoQyKPMCore.kpm_mul — Function
kpm_mul(
+) where {T<:AbstractFloat}
Given the KPM moments $\mu$, calculate the spectral density $\rho(\epsilon)$ and corresponding energies $\epsilon$ where it is evaluated. Here $W$ is the spectral bandwidth.
source
SmoQyKPMCore.kpm_dot — Function
kpm_dot(
+    A, coefs::AbstractVector, bounds, R::T,
+    tmp = zeros(eltype(v), size(v)..., 3)
+) where {T<:AbstractVecOrMat}
If $R$ is a single vector, then calculate the inner product
\[\begin{align*}
+S & = \langle R | F(A) | R \rangle \\
+S & = \sum_{m=1}^M \langle R | c_m T_m(A^\prime) | R \rangle
+\end{align*},\]
wher $A^\prime$ is the scaled version of $A$. If $R$ is a matrix, then calculate
\[\begin{align*}
+S & = \langle R | F(A) | R \rangle \\
+S & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle R_n | c_m T_m(A^\prime) | R_n \rangle
+\end{align*},\]
where $| R_n \rangle$ is a column of $R$.
source
kpm_dot(
+    A, coefs::AbstractVector, bounds, U::T, V::T,
+    tmp = zeros(eltype(v), size(v)..., 3)
+) where {T<:AbstractVecOrMat}
If $U$ and $V$ are single vectors, then calculate the inner product
\[\begin{align*}
+S & = \langle U | F(A) | V \rangle \\
+  & = \sum_{m=1}^M \langle U | c_m T_m(A^\prime) | V \rangle
+\end{align*},\]
wher $A^\prime$ is the scaled version of $A$. If $U$ and $V$ are matrices, then calculate
\[\begin{align*}
+S & = \langle U | F(A) | V \rangle \\
+  & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle U_n | c_m T_m(A^\prime) | V_n \rangle
+\end{align*},\]
where $| U_n \rangle$ and $| V_n \rangle$ are the columns of each matrix.
source
kpm_dot(
+    A, kpm_expansion::KPMExpansion, R::T,
+    tmp = zeros(eltype(R), size(R)..., 3)
+) where {T<:AbstractVecOrMat}
If $R$ is a single vector, then calculate the inner product
\[\begin{align*}
+S & = \langle R | F(A) | R \rangle \\
+S & = \sum_{m=1}^M \langle R | c_m T_m(A^\prime) | R \rangle
+\end{align*},\]
wher $A^\prime$ is the scaled version of $A$. If $R$ is a matrix, then calculate
\[\begin{align*}
+S & = \langle R | F(A) | R \rangle \\
+S & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle R_n | c_m T_m(A^\prime) | R_n \rangle
+\end{align*},\]
where $| R_n \rangle$ is a column of $R$.
source
kpm_dot(
+    A, kpm_expansion::KPMExpansion, U::T, V::T,
+    tmp = zeros(eltype(V), size(V)..., 3)
+) where {T<:AbstractVecOrMat}
If $U$ and $V$ are single vectors, then calculate the inner product
\[\begin{align*}
+S & = \langle U | F(A) | V \rangle \\
+  & = \sum_{m=1}^M \langle U | c_m T_m(A^\prime) | V \rangle
+\end{align*},\]
wher $A^\prime$ is the scaled version of $A$. If $U$ and $V$ are matrices, then calculate
\[\begin{align*}
+S & = \langle U | F(A) | V \rangle \\
+  & = \frac{1}{N} \sum_{n=1}^N \sum_{m=1}^M \langle U_n | c_m T_m(A^\prime) | V_n \rangle
+\end{align*},\]
where $| U_n \rangle$ and $| V_n \rangle$ are the columns of each matrix.
source
SmoQyKPMCore.kpm_mul — Function
kpm_mul(
     A, coefs::AbstractVector, bounds, v::T, tmp = zeros(eltype(v), size(v)..., 3)
-) where {T<:AbstractVecOrMat}
Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.
source
kpm_mul(A, kpm_expansion::KPMExpansion, v::T) where {T<:AbstractVector}
Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.
source
SmoQyKPMCore.kpm_mul! — Function
kpm_mul!(
+) where {T<:AbstractVecOrMat}
Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.
source
kpm_mul(A, kpm_expansion::KPMExpansion, v::T) where {T<:AbstractVecOrMat}
Evaluate and return the vector $v^\prime = F(A) \cdot v$ where $F(A)$ is represented by the Chebyshev expansion. For more information refer to kpm_mul!.
source
SmoQyKPMCore.kpm_mul! — Function
kpm_mul!(
     v′::T, A, coefs::AbstractVector, bounds, v::T,
     tmp = zeros(eltype(v), size(v)..., 3)
-) where {T<:AbstractVecOrMat}
Evaluates $v^\prime = F(A) \cdot v$, writing the result to v′, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate $F(A)$, where the eigenspectrum of $A$ is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for $F(A)$. Lastly, tmp is an array used to avoid dynamic memory allocations.
source
kpm_mul!(
+) where {T<:AbstractVecOrMat}
Evaluates $v^\prime = F(A) \cdot v$, writing the result to v′, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate $F(A)$, where the eigenspectrum of $A$ is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for $F(A)$. Lastly, tmp is an array used to avoid dynamic memory allocations.
source
kpm_mul!(
     v′::T, A, kpm_expansion::KPMExpansion, v::T,
     tmp = zeros(eltype(v), size(v)..., 3)
-) where {T<:AbstractVecOrMat}
Evaluates $v^\prime = F(A) \cdot v$, writing the result to v′, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.
source
SmoQyKPMCore.kpm_eval — Function
kpm_eval(x::AbstractFloat, coefs, bounds)
Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(A::AbstractMatrix, coefs, bounds)
Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(x::T, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}
Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(A::AbstractMatrix{T}, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}
Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
SmoQyKPMCore.kpm_eval! — Function
kpm_eval!(
+) where {T<:AbstractVecOrMat}
Evaluates $v^\prime = F(A) \cdot v$, writing the result to v′, where $F(A)$ is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate $u = A\cdot v$, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.
source
SmoQyKPMCore.kpm_eval — Function
kpm_eval(x::AbstractFloat, coefs, bounds)
Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(A::AbstractMatrix, coefs, bounds)
Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(x::T, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}
Evaluate $F(x)$ where $x$ is real number in the interval bounds[1] < x < bound[2], and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
kpm_eval(A::AbstractMatrix{T}, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}
Evaluate and return the matrix $F(A),$ where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs.
source
SmoQyKPMCore.kpm_eval! — Function
kpm_eval!(
     F::AbstractMatrix, A, coefs::AbstractVector, bounds,
     tmp = zeros(eltype(F), size(F)..., 3)
-)
Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, tmp is used to avoid dynamic memory allocations.
source
kpm_eval!(
+)
Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, tmp is used to avoid dynamic memory allocations.
source
kpm_eval!(
     F::AbstractMatrix, A, kpm_expansion::KPMExpansion,
     tmp = zeros(eltype(F), size(F)..., 3)
-)
Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, the array tmp is used to avoid dynamic memory allocations.
source
SmoQyKPMCore.apply_jackson_kernel! — Function
apply_jackson_kernel!(coefs)
Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.
source
apply_jackson_kernel!(kpm_expansion::KPMExpansion)
Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.
source
SmoQyKPMCore.lanczos — Function
lanczos(niters, v, A, S = I, rng = Random.default_rng())
Use niters Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.
This function returns a SymTridiagonal matrix. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.
Similar generalizations of Lanczos have been considered in [2] and [3].
1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
https://doi.org/10.1093%2Fimamat%2F10.3.373.
2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
https://doi.org/10.1090/s0025-5718-1982-0669648-0
3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
https://doi.org/10.1016/j.commatsci.2011.02.021.
source
SmoQyKPMCore.lanczos! — Function
lanczos!(
+)
Evaluate and write the matrix $F(A)$ to F, where $A$ is an operator with strictly real eigenvalues and the function $F(\bullet)$ is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, the array tmp is used to avoid dynamic memory allocations.
source
SmoQyKPMCore.apply_jackson_kernel! — Function
apply_jackson_kernel!(coefs)
Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.
source
apply_jackson_kernel!(kpm_expansion::KPMExpansion)
Modify the Chebyshev expansion coefficients by applying the Jackson kernel to them.
source
SmoQyKPMCore.lanczos — Function
lanczos(niters, v, A, S = I, rng = Random.default_rng())
Use niters Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.
This function returns a SymTridiagonal matrix. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.
Similar generalizations of Lanczos have been considered in [2] and [3].
1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
https://doi.org/10.1093%2Fimamat%2F10.3.373.
2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
https://doi.org/10.1090/s0025-5718-1982-0669648-0
3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
https://doi.org/10.1016/j.commatsci.2011.02.021.
source
SmoQyKPMCore.lanczos! — Function
lanczos!(
     αs::AbstractVector, βs::AbstractVector, v::AbstractVector,
     A, S = I,
     tmp::AbstractMatrix = zeros(eltype(v), length(v), 5);
@@ -57,4 +97,4 @@
     A, S = I,
     tmp::AbstractArray = zeros(eltype(v), size(v)..., 5);
     rng = Random.default_rng()
-)
Use Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.
The number of Lanczos iterations performed equals niters = length(αs), and niters - 1 == length(βs). This function returns a SymTridiagonal matrix based on the contents of the vectors αs and βs. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.
Note that if αs, βs and v are all matrices, then each column of v is treated as a seperate vector, and a vector of SymTridiagonal of length size(v,2) will be returned.
Similar generalizations of Lanczos have been considered in [2] and [3].
1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
https://doi.org/10.1093%2Fimamat%2F10.3.373.
2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
https://doi.org/10.1090/s0025-5718-1982-0669648-0
3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
https://doi.org/10.1016/j.commatsci.2011.02.021.
source

+)

Use Lanczos iterations to find a truncated tridiagonal representation of $A\cdot S$, up to similarity transformation. Here, $A$ is any Hermitian matrix, while $S$ is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, $S = I$. The extension to non-identity matrices $S$ is as follows: Each matrix-vector product $A\cdot v$ becomes $(A S) \cdot v$, and each vector inner product $w^\dagger \cdot v$ becomes $w^\dagger \cdot S \cdot v$. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for $A$ and $S$, respectively.

The number of Lanczos iterations performed equals niters = length(αs), and niters - 1 == length(βs). This function returns a SymTridiagonal matrix based on the contents of the vectors αs and βs. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.

Note that if αs, βs and v are all matrices, then each column of v is treated as a seperate vector, and a vector of SymTridiagonal of length size(v,2) will be returned.

Similar generalizations of Lanczos have been considered in [2] and [3].

1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),

https://doi.org/10.1093%2Fimamat%2F10.3.373.

2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),

https://doi.org/10.1090/s0025-5718-1982-0669648-0

3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),

https://doi.org/10.1016/j.commatsci.2011.02.021.

source
diff --git a/dev/index.html b/dev/index.html index 8430f34..7161aff 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · SmoQyKPMCore.jl
GitHub

SmoQyKPMCore

Documentation for SmoQyKPMCore package. The SmoQyKPMCore package implements and exports an optimized, low-level implementation of the Kernel Polynomial Method (KPM) algorithm for approximating functions of operators with strictly real, bounded eigenvalues via a Chebyshev polynomial expansion.

Funding

The development of this code was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0022311.

+Home · SmoQyKPMCore.jl
GitHub

SmoQyKPMCore

Documentation for SmoQyKPMCore package. The SmoQyKPMCore package implements and exports an optimized, low-level implementation of the Kernel Polynomial Method (KPM) algorithm for approximating functions of operators with strictly real, bounded eigenvalues via a Chebyshev polynomial expansion.

Funding

The development of this code was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0022311.

diff --git a/dev/objects.inv b/dev/objects.inv index db3ea828293f3cb2a79b6b600b259bed8e689630..1488468a84a72eb0b460aafec9914754957ec033 100644 GIT binary patch delta 410 zcmV;L0cHM@1eXMmhksLDZ-OusedkwnnY_uIW|@g^PG@Eo6Bp+*skAs0dK)MUY}tQb zfoebj%Tr4^=iGC@JVt1OB0(ku9Apgm^pa_mBCGxIT}i{71)sncFc*alDz|&Uvyjfiqb9%A)#R>)tu2F>Hu}4Ir3fv&2_7kSd9BD3O+|<0RuTAm(0TOdQ EDRYX#WB>pF delta 405 zcmV;G0c!r21d;@hhkuh%Z-Ouo$KUxBT_$fbr&(s=o70(@#l*$=Oe!r7gg0$cXo zSD+eDKzM2?_y67hu04+t8ly;%2>}OLLxkQSpoWa_o)(X@g8|y6i-EaU{HHz#Xu4h+ zn*ctN_w{N9VJ`OZj{0fHeGgg2oU+Ihi3P5;5m?a3nsOdO?0*GVc)-evJfYHQ6#^IK zDpXP#Dsfqv+R5^5&qZ8QEn^`W@z}n>=L(C^i#wUjB^M5hZF#_)`-y$P&M^XJA8E*_ ze}@kK#Ffm+GUy6!s?5<)nJtAKd0sQRMYq+NJU{8>FJ%T)d;zI6 zc}+9NFDITknwg_V#mqMuu(%!1FU|RwIX}>$nX-G)sX!#0HiC&Bg|^-#u0iQp=<+sN4i1^DyO8v*>hk>q1hX%dbK}OM diff --git a/dev/search_index.js b/dev/search_index.js index 1d8284a..071ca8d 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"api/#API","page":"API","title":"API","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"Note that for many methods there are two verions, one that relies on taking an instance of the KPMExpansion type as an argument, and a lower level one that does not.","category":"page"},{"location":"api/","page":"API","title":"API","text":"KPMExpansion\nkpm_update!\nkpm_update_bounds!\nkpm_update_order!\nkpm_coefs\nkpm_coefs!\nkpm_moments\nkpm_moments!\nkpm_density\nkpm_density!\nkpm_mul\nkpm_mul!\nkpm_eval\nkpm_eval!\napply_jackson_kernel\napply_jackson_kernel!\nlanczos\nlanczos!","category":"page"},{"location":"api/","page":"API","title":"API","text":"KPMExpansion\nKPMExpansion(::Function, ::Any, ::Int, ::Int)\nkpm_update!\nkpm_update_bounds!\nkpm_update_order!\nkpm_coefs\nkpm_coefs!\nkpm_moments\nkpm_moments!\nkpm_density\nkpm_density!\nkpm_mul\nkpm_mul!\nkpm_eval\nkpm_eval!\napply_jackson_kernel\napply_jackson_kernel!\nlanczos\nlanczos!","category":"page"},{"location":"api/#SmoQyKPMCore.KPMExpansion","page":"API","title":"SmoQyKPMCore.KPMExpansion","text":"mutable struct KPMExpansion{T<:AbstractFloat, Tfft<:FFTW.r2rFFTWPlan}\n\nA type to represent the Chebyshev polynomial expansion used in the KPM algorithm.\n\nFields\n\nM::Int: Order of the Chebyshev polynomial expansion.\nbounds::NTuple{2,T}: Bounds on eigenspectrum in the KPM algorithm.\nbuf::Vector{T}: The first M elements of this vector of the Chebyshev expansion coefficients.\nr2rplan::Tfft: Plan for performing DCT to efficiently calculate the Chebyshev expansion coefficients via Chebyshev-Gauss quadrature.\n\n\n\n\n\n","category":"type"},{"location":"api/#SmoQyKPMCore.KPMExpansion-Tuple{Function, Any, Int64, Int64}","page":"API","title":"SmoQyKPMCore.KPMExpansion","text":"KPMExpansion(func::Function, bounds, M::Int, N::Int = 2*M)\n\nInitialize an instance of the KPMExpansion type to approximate the univariate function func, called as func(x), with a order M Chebyshev polynomial expansion on the interval bounds[1] < x bounds[2]. Here, N ≥ M is the number of points at which func is evaluated on that specified interval, which are then used to calculate the expansion coeffiencents via Chebyshev-Gauss quadrature.\n\n\n\n\n\n","category":"method"},{"location":"api/#SmoQyKPMCore.kpm_update!","page":"API","title":"SmoQyKPMCore.kpm_update!","text":"kpm_update!(\n kpm_expansion::KPMExpansion{T}, func::Function, bounds, M::Int, N::Int = 2*M\n) where {T<:AbstractFloat}\n\nIn-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, expansion order M, and number of points at which the expanded function is evaluated when computing the expansion coefficients. This includes recomputing the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_update_bounds!","page":"API","title":"SmoQyKPMCore.kpm_update_bounds!","text":"kpm_update_bounds!(\n kpm_expansion::KPMExpansion{T}, func::Function, bounds\n) where {T<:AbstractFloat}\n\nIn-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, recomputing the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_update_order!","page":"API","title":"SmoQyKPMCore.kpm_update_order!","text":"kpm_update_order!(\n kpm_expansion::KPMExpansion, func::Function, M::Int, N::Int = 2*M\n)\n\nIn-place update the expansion order M for an instance of KPMExpansion, recomputing the expansion coefficients. It is also possible to udpate the number of point N the function func is evaluated at to calculate the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_coefs","page":"API","title":"SmoQyKPMCore.kpm_coefs","text":"kpm_coefs(func::Function, bounds, M::Int, N::Int = 2*M)\n\nCalculate and return the Chebyshev expansion coefficients. Refer to kpm_coefs! for more information.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_coefs!","page":"API","title":"SmoQyKPMCore.kpm_coefs!","text":"kpm_coefs!(\n coefs::AbstractVector{T}, func::Function, bounds,\n buf::Vector{T} = zeros(T, 2*length(coefs)),\n r2rplan = FFTW.plan_r2r!(zeros(T, 2*length(coefs)), FFTW.REDFT10)\n) where {T<:AbstractFloat}\n\nCalculate and record the Chebyshev polynomial expansion coefficients to order M in the vector ceofs for the function func on the interval (bounds[1], bounds[2]). Let length(buf) be the number of evenly spaced points on the interval for which func is evaluated when performing Chebyshev-Gauss quadrature to compute the Chebyshev polynomial expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_moments","page":"API","title":"SmoQyKPMCore.kpm_moments","text":"kpm_moments(\n M::Int, A, bounds, R::T, tmp = zeros(eltype(R), size(R)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate and return the first M moments\n\nmu_m = langle R T_m(A) R rangle\n\nif R is a vector, and if R is a matrix then\n\nmu_m = frac1N sum_n=1^N langle R_n T_m(A) R_n rangle\n\nwhere R_n rangle is the n'th column of R.\n\n\n\n\n\nkpm_moments(\n M::Int, A, bounds, U::T, V::T, tmp = zeros(eltype(V), size(V)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate and return the first M moments\n\nmu_m = langle U T_m(A) V rangle\n\nif U and V are vector, and if they are matrices then\n\nmu_m = frac1N sum_n=1^N langle U_n T_m(A) V_n rangle\n\nwhere U_n rangle and V_n rangle are are the n'th columns of each matrix.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_moments!","page":"API","title":"SmoQyKPMCore.kpm_moments!","text":"kpm_moments!(\n μ::AbstractVector, A, bounds, R::T, tmp = zeros(eltype(R), size(R)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate the moments\n\nmu_m = langle R T_m(A) R rangle\n\nif R is a vector, and if R is a matrix then\n\nmu_m = frac1N sum_n=1^N langle R_n T_m(A) R_n rangle\n\nwhere R_n rangle is the n'th column of R.\n\n\n\n\n\nkpm_moments!(\n μ::AbstractVector, A, bounds, U::T, V::T, tmp = zeros(eltype(V), size(V)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate the moments\n\nmu_m = langle U T_m(A) V rangle\n\nif U and V are vector, and if they are matrices then\n\nmu_m = frac1N sum_n=1^N langle U_n T_m(A) V_n rangle\n\nwhere U_n rangle and V_n rangle are are the n'th columns of each matrix.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_density","page":"API","title":"SmoQyKPMCore.kpm_density","text":"kpm_density(\n N::Int, μ::T, W\n) where {T<:AbstractVector}\n\nGiven the KPM moments mu, evaluate the corresponding spectral density rho(epsilon) at N energy points epsilon, where W is the spectral bandwidth or the bounds of the spectrum. This function returns both rho(epsilon) and epsilon.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_density!","page":"API","title":"SmoQyKPMCore.kpm_density!","text":"kpm_density!(\n ρ::AbstractVector{T},\n ϵ::AbstractVector{T},\n μ::AbstractVector{T},\n bounds,\n r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01)\n) where {T<:AbstractFloat}\n\nGiven the KPM moments mu, calculate the spectral density rho(epsilon) and corresponding energies epsilon where it is evaluated. Here bounds is bounds the eigenspecturm, such that the bandwidth is given by W = bounds[2] - bounds[1].\n\n\n\n\n\nkpm_density!(\n ρ::AbstractVector{T},\n ϵ::AbstractVector{T},\n μ::AbstractVector{T},\n W::T,\n r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01),\n) where {T<:AbstractFloat}\n\nGiven the KPM moments mu, calculate the spectral density rho(epsilon) and corresponding energies epsilon where it is evaluated. Here W is the spectral bandwidth.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_mul","page":"API","title":"SmoQyKPMCore.kpm_mul","text":"kpm_mul(\n A, coefs::AbstractVector, bounds, v::T, tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluate and return the vector v^prime = F(A) cdot v where F(A) is represented by the Chebyshev expansion. For more information refer to kpm_mul!.\n\n\n\n\n\nkpm_mul(A, kpm_expansion::KPMExpansion, v::T) where {T<:AbstractVector}\n\nEvaluate and return the vector v^prime = F(A) cdot v where F(A) is represented by the Chebyshev expansion. For more information refer to kpm_mul!.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_mul!","page":"API","title":"SmoQyKPMCore.kpm_mul!","text":"kpm_mul!(\n v′::T, A, coefs::AbstractVector, bounds, v::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluates v^prime = F(A) cdot v, writing the result to v′, where F(A) is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate u = Acdot v, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate F(A), where the eigenspectrum of A is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for F(A). Lastly, tmp is an array used to avoid dynamic memory allocations.\n\n\n\n\n\nkpm_mul!(\n v′::T, A, kpm_expansion::KPMExpansion, v::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluates v^prime = F(A) cdot v, writing the result to v′, where F(A) is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate u = Acdot v, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_eval","page":"API","title":"SmoQyKPMCore.kpm_eval","text":"kpm_eval(x::AbstractFloat, coefs, bounds)\n\nEvaluate F(x) where x is real number in the interval bounds[1] < x < bound[2], and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(A::AbstractMatrix, coefs, bounds)\n\nEvaluate and return the matrix F(A) where A is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(x::T, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}\n\nEvaluate F(x) where x is real number in the interval bounds[1] < x < bound[2], and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(A::AbstractMatrix{T}, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}\n\nEvaluate and return the matrix F(A) where A is an operator with strictly real eigenvalues and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_eval!","page":"API","title":"SmoQyKPMCore.kpm_eval!","text":"kpm_eval!(\n F::AbstractMatrix, A, coefs::AbstractVector, bounds,\n tmp = zeros(eltype(F), size(F)..., 3)\n)\n\nEvaluate and write the matrix F(A) to F, where A is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\nkpm_eval!(\n F::AbstractMatrix, A, kpm_expansion::KPMExpansion,\n tmp = zeros(eltype(F), size(F)..., 3)\n)\n\nEvaluate and write the matrix F(A) to F, where A is an operator with strictly real eigenvalues and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, the array tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.apply_jackson_kernel","page":"API","title":"SmoQyKPMCore.apply_jackson_kernel","text":"apply_jackson_kernel(coefs)\n\nReturn the Chebyshev expansion coefficients transformed by the Jackson kernel.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.apply_jackson_kernel!","page":"API","title":"SmoQyKPMCore.apply_jackson_kernel!","text":"apply_jackson_kernel!(coefs)\n\nModify the Chebyshev expansion coefficients by applying the Jackson kernel to them.\n\n\n\n\n\napply_jackson_kernel!(kpm_expansion::KPMExpansion)\n\nModify the Chebyshev expansion coefficients by applying the Jackson kernel to them.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.lanczos","page":"API","title":"SmoQyKPMCore.lanczos","text":"lanczos(niters, v, A, S = I, rng = Random.default_rng())\n\nUse niters Lanczos iterations to find a truncated tridiagonal representation of Acdot S, up to similarity transformation. Here, A is any Hermitian matrix, while S is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, S = I. The extension to non-identity matrices S is as follows: Each matrix-vector product Acdot v becomes (A S) cdot v, and each vector inner product w^dagger cdot v becomes w^dagger cdot S cdot v. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for A and S, respectively.\n\nThis function returns a SymTridiagonal matrix. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.\n\nSimilar generalizations of Lanczos have been considered in [2] and [3].\n\n\n\n1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),\n\nhttps://doi.org/10.1093%2Fimamat%2F10.3.373.\n\n2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),\n\nhttps://doi.org/10.1090/s0025-5718-1982-0669648-0\n\n3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),\n\nhttps://doi.org/10.1016/j.commatsci.2011.02.021.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.lanczos!","page":"API","title":"SmoQyKPMCore.lanczos!","text":"lanczos!(\n αs::AbstractVector, βs::AbstractVector, v::AbstractVector,\n A, S = I,\n tmp::AbstractMatrix = zeros(eltype(v), length(v), 5);\n rng = Random.default_rng()\n)\n\nlanczos!(\n αs::AbstractMatrix, βs::AbstractMatrix, v::AbstractMatrix,\n A, S = I,\n tmp::AbstractArray = zeros(eltype(v), size(v)..., 5);\n rng = Random.default_rng()\n)\n\nUse Lanczos iterations to find a truncated tridiagonal representation of Acdot S, up to similarity transformation. Here, A is any Hermitian matrix, while S is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, S = I. The extension to non-identity matrices S is as follows: Each matrix-vector product Acdot v becomes (A S) cdot v, and each vector inner product w^dagger cdot v becomes w^dagger cdot S cdot v. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for A and S, respectively.\n\nThe number of Lanczos iterations performed equals niters = length(αs), and niters - 1 == length(βs). This function returns a SymTridiagonal matrix based on the contents of the vectors αs and βs. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.\n\nNote that if αs, βs and v are all matrices, then each column of v is treated as a seperate vector, and a vector of SymTridiagonal of length size(v,2) will be returned.\n\nSimilar generalizations of Lanczos have been considered in [2] and [3].\n\n\n\n1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),\n\nhttps://doi.org/10.1093%2Fimamat%2F10.3.373.\n\n2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),\n\nhttps://doi.org/10.1090/s0025-5718-1982-0669648-0\n\n3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),\n\nhttps://doi.org/10.1016/j.commatsci.2011.02.021.\n\n\n\n\n\n","category":"function"},{"location":"usage/","page":"Usage","title":"Usage","text":"EditURL = \"../../examples/usage.jl\"","category":"page"},{"location":"usage/#Usage","page":"Usage","title":"Usage","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"Here we demonstrate the basic usage for the SmoQyKPMCore package. Let us first import the relevant packages we will want to use in this example.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"using LinearAlgebra\nusing SparseArrays","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Package for making figures","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"using CairoMakie\nCairoMakie.activate!(type = \"svg\")\n\nusing SmoQyKPMCore","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"In this usage example we will consider a 1D chain tight-binding model","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"hatH = -t sum_isigma (hatc_i+1sigma^dagger hatc_isigma + rm Hc)\n = sum_sigma hatc_isigma^dagger H_ij hatc_jsigma","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"where hatc_isigma^dagger (hatc_isigma) creates (annihilates) a spin-sigma electron on site i in the lattice, and t is the nearest-neighbor hopping amplitude.","category":"page"},{"location":"usage/#Density-Matrix-Approximation","page":"Usage","title":"Density Matrix Approximation","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"With H the Hamiltonian matrix, the corresponding density matrix is given by rho = f(H) where","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"f(epsilon) = frac11+e^beta (epsilon-mu) = frac12 left 1 + tanhleft(tfracbeta(epsilon-mu)2right) right","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"is the Fermi function, with beta = 1T the inverse temperature and mu the chemical potential. Here we will applying the kernel polynomial method (KPM) algorithm to approximate the density matrix rho by a Chebyshev polynomial expansion.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us first define the relevant model parameter values.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Nearest-neighbor hopping amplitude.\nt = 1.0\n\n# Chemical potential.\nμ = 0.0\n\n# System size.\nL = 16\n\n# Inverse temperature.\nβ = 4.0;\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now we need to construct the Hamiltonian matrix H.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"H = zeros(L,L)\nfor i in 1:L\n j = mod1(i+1,L)\n H[j,i], H[i,j] = -t, -t\nend\nH","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"We also need to define the Fermi function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Fermi function.\nfermi(ϵ, μ, β) = (1+tanh(β*(ϵ-μ)/2))/2","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us calculate the exact density matrix rho so that we may asses the accuracy of the KPM method.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Diagonalize the Hamiltonian matrix.\nϵ, U = eigen(H)\n\n# Calculate the density matrix.\nρ = U * Diagonal(fermi.(ϵ, μ, β)) * adjoint(U)","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now let use the KPMExpansion type to approximate rho. We will need to define the order M of the expansion and give approximate bounds for the eigenspectrum of H, making sure to overestimate the the true interval spanned by the eigenvalues of H. Because we are considering the simple non-interacting model, the exact eigenspectrum of H is known and spans the interval epsilon in -2t 2t","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Define eigenspectrum bounds.\nbounds = (-3.0t, 3.0t)\n\n# Define order of Chebyshev expansion used in KPM approximation.\nM = 10\n\n# Initialize expansion.\nkpm_expansion = KPMExpansion(x -> fermi(x, μ, β), bounds, M);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now let us test and see how good a job it does approximating the density matrix rho, using the kpm_eval! function to do so.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Initialize KPM density matrix approxiatmion\nρ_kpm = similar(H)\n\n# Initialize additional matrices to avoid dynamic memory allocation.\nmtmp = zeros(L,L,3)\n\n# Calculate KPM density matrix approximation.\nkpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)\n\n# Check how good an approximation it is.\nprintln(\"Matrix Error = \", norm(ρ_kpm - ρ)/norm(ρ) )","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"It is also possible to efficiently multiply vectors by KPM approximation to the density matrix rho. Let us test this functionality with the kpm_mul! on a random vector, seeing how accurate the result is.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Initialize random vector.\nv = randn(L)\n\n# Calculate exact product with density matrix.\nρv = ρ * v\n\n# Initialize vector to contain approximate product with densith matrix.\nρv_kpm = similar(v)\n\n# Initialize to avoid dynamic memory allocation.\nvtmp = zeros(L, 3)\n\n# Calculate approximate product with density matrix.\nkpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)\n\n# Check how good the approximation is.\nprintln(\"Vector Error = \", norm(ρv_kpm - ρv) / norm(ρv) )","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us now provide the KPM approximation with better bounds on the eigenspectrum and also increase the order the of the expansion and see how the result improves. We will use the kpm_update! function to update the KPMExpansion in-place.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Define eigenspectrum bounds.\nbounds = (-2.5t, 2.5t)\n\n# Define order of Chebyshev expansion used in KPM approximation.\nM = 100\n\n# Update KPM approximation.\nkpm_update!(kpm_expansion, x -> fermi(x, μ, β), bounds, M)\n\n# Calculate KPM density matrix approximation.\nkpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)\n\n# Check how good an approximation it is.\nprintln(\"Matrix Error = \", norm(ρ_kpm - ρ)/norm(ρ) )\n\n# Calculate approximate product with density matrix.\nkpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)\n\n# Check how good the approximation is.\nprintln(\"Vector Error = \", norm(ρv_kpm - ρv) / norm(ρv) )","category":"page"},{"location":"usage/#Density-of-States-Approximation","page":"Usage","title":"Density of States Approximation","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"Next let us demonstrate how we may approximate the density of states mathcalN(epsilon) for a 1D chain tight-binding model. The first step a very larger Hamiltonian matrix H which we represent as a spare matrix.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Size of 1D chain considered\nL = 10_000\n\n# Construct sparse Hamiltonian matrix.\nrows = Int[]\ncols = Int[]\nfor i in 1:L\n j = mod1(i+1, L)\n push!(rows,i)\n push!(cols,j)\n push!(rows,j)\n push!(cols,i)\nend\nvals = fill(-t, length(rows))\nH = sparse(rows, cols, vals)","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us now calculate the fist M moments mu_m using N random vectors using the function kpm_moments.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Number of moments to calculate.\nM = 250\n\n# Number of random vectors used to approximate moments.\nN = 100\n\n# Initialize random vectors.\nR = randn(L, N)\n\n# Calculate the moments.\nμ = kpm_moments(M, H, bounds, R);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Having calculate the moments, let us next evaluate the density of states mathcalN(epsilon) at P points with the kpm_density function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Number of points at which to evaluate the density of states.\nP = 1000\n\n# Evaluate density of states.\ndos, ϵ = kpm_density(P, μ, bounds);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Without regularization, the approximate for the density of states generated above will have clearly visible Gibbs oscillations. To partially suppress these artifacts, we apply the Jackson kernel to the moments mu using the apply_jackson_kernel function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Apply Jackson kernel.\nμ_jackson = apply_jackson_kernel(μ)\n\n# Evaluate density of states.\ndos_jackson, ϵ_jackson = kpm_density(P, μ_jackson, bounds);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Have approximated the density of states, let us now plot the result.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"fig = Figure(\n size = (700, 400),\n fonts = (; regular= \"CMU Serif\"),\n figure_padding = 10\n)\n\nax = Axis(\n fig[1, 1],\n aspect = 7/4,\n xlabel = L\"\\epsilon\", ylabel = L\"\\mathcal{N}(\\epsilon)\",\n xticks = range(start = bounds[1], stop = bounds[2], length = 5),\n xlabelsize = 30, ylabelsize = 30,\n xticklabelsize = 24, yticklabelsize = 24,\n)\n\nlines!(\n ϵ, dos,\n linewidth = 2.0,\n alpha = 1.0,\n color = :red,\n linestyle = :solid,\n label = \"Dirichlet\"\n)\n\nlines!(\n ϵ_jackson, dos_jackson,\n linewidth = 3.0,\n alpha = 1.0,\n color = :black,\n linestyle = :solid,\n label = \"Jackson\"\n)\n\nxlims!(\n ax, bounds[1], bounds[2]\n)\n\nylims!(\n ax, 0.0, 1.05 * maximum(dos)\n)\n\naxislegend(\n ax, halign = :center, valign = :top, labelsize = 30\n)\n\nfig","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = SmoQyKPMCore","category":"page"},{"location":"#SmoQyKPMCore","page":"Home","title":"SmoQyKPMCore","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Documentation for SmoQyKPMCore package. The SmoQyKPMCore package implements and exports an optimized, low-level implementation of the Kernel Polynomial Method (KPM) algorithm for approximating functions of operators with strictly real, bounded eigenvalues via a Chebyshev polynomial expansion.","category":"page"},{"location":"#Funding","page":"Home","title":"Funding","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The development of this code was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0022311.","category":"page"}] +[{"location":"api/#API","page":"API","title":"API","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"Note that for many methods there are two verions, one that relies on taking an instance of the KPMExpansion type as an argument, and a lower level one that does not.","category":"page"},{"location":"api/","page":"API","title":"API","text":"KPMExpansion\nkpm_update!\nkpm_update_bounds!\nkpm_update_order!\nkpm_coefs\nkpm_coefs!\nkpm_moments\nkpm_moments!\nkpm_density\nkpm_density!\nkpm_dot\nkpm_mul\nkpm_mul!\nkpm_eval\nkpm_eval!\napply_jackson_kernel\napply_jackson_kernel!\nlanczos\nlanczos!","category":"page"},{"location":"api/","page":"API","title":"API","text":"KPMExpansion\nKPMExpansion(::Function, ::Any, ::Int, ::Int)\nkpm_update!\nkpm_update_bounds!\nkpm_update_order!\nkpm_coefs\nkpm_coefs!\nkpm_moments\nkpm_moments!\nkpm_density\nkpm_density!\nkpm_dot\nkpm_mul\nkpm_mul!\nkpm_eval\nkpm_eval!\napply_jackson_kernel\napply_jackson_kernel!\nlanczos\nlanczos!","category":"page"},{"location":"api/#SmoQyKPMCore.KPMExpansion","page":"API","title":"SmoQyKPMCore.KPMExpansion","text":"mutable struct KPMExpansion{T<:AbstractFloat, Tfft<:FFTW.r2rFFTWPlan}\n\nA type to represent the Chebyshev polynomial expansion used in the KPM algorithm.\n\nFields\n\nM::Int: Order of the Chebyshev polynomial expansion.\nbounds::NTuple{2,T}: Bounds on eigenspectrum in the KPM algorithm.\nbuf::Vector{T}: The first M elements of this vector of the Chebyshev expansion coefficients.\nr2rplan::Tfft: Plan for performing DCT to efficiently calculate the Chebyshev expansion coefficients via Chebyshev-Gauss quadrature.\n\n\n\n\n\n","category":"type"},{"location":"api/#SmoQyKPMCore.KPMExpansion-Tuple{Function, Any, Int64, Int64}","page":"API","title":"SmoQyKPMCore.KPMExpansion","text":"KPMExpansion(func::Function, bounds, M::Int, N::Int = 2*M)\n\nInitialize an instance of the KPMExpansion type to approximate the univariate function func, called as func(x), with a order M Chebyshev polynomial expansion on the interval bounds[1] < x bounds[2]. Here, N ≥ M is the number of points at which func is evaluated on that specified interval, which are then used to calculate the expansion coeffiencents via Chebyshev-Gauss quadrature.\n\n\n\n\n\n","category":"method"},{"location":"api/#SmoQyKPMCore.kpm_update!","page":"API","title":"SmoQyKPMCore.kpm_update!","text":"kpm_update!(\n kpm_expansion::KPMExpansion{T}, func::Function, bounds, M::Int, N::Int = 2*M\n) where {T<:AbstractFloat}\n\nIn-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, expansion order M, and number of points at which the expanded function is evaluated when computing the expansion coefficients. This includes recomputing the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_update_bounds!","page":"API","title":"SmoQyKPMCore.kpm_update_bounds!","text":"kpm_update_bounds!(\n kpm_expansion::KPMExpansion{T}, func::Function, bounds\n) where {T<:AbstractFloat}\n\nIn-place update an instance of KPMExpansion to reflect new values for eigenspectrum bounds, recomputing the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_update_order!","page":"API","title":"SmoQyKPMCore.kpm_update_order!","text":"kpm_update_order!(\n kpm_expansion::KPMExpansion, func::Function, M::Int, N::Int = 2*M\n)\n\nIn-place update the expansion order M for an instance of KPMExpansion, recomputing the expansion coefficients. It is also possible to udpate the number of point N the function func is evaluated at to calculate the expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_coefs","page":"API","title":"SmoQyKPMCore.kpm_coefs","text":"kpm_coefs(func::Function, bounds, M::Int, N::Int = 2*M)\n\nCalculate and return the Chebyshev expansion coefficients. Refer to kpm_coefs! for more information.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_coefs!","page":"API","title":"SmoQyKPMCore.kpm_coefs!","text":"kpm_coefs!(\n coefs::AbstractVector{T}, func::Function, bounds,\n buf::Vector{T} = zeros(T, 2*length(coefs)),\n r2rplan = FFTW.plan_r2r!(zeros(T, 2*length(coefs)), FFTW.REDFT10)\n) where {T<:AbstractFloat}\n\nCalculate and record the Chebyshev polynomial expansion coefficients to order M in the vector ceofs for the function func on the interval (bounds[1], bounds[2]). Let length(buf) be the number of evenly spaced points on the interval for which func is evaluated when performing Chebyshev-Gauss quadrature to compute the Chebyshev polynomial expansion coefficients.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_moments","page":"API","title":"SmoQyKPMCore.kpm_moments","text":"kpm_moments(\n M::Int, A, bounds, R::T,\n tmp = zeros(eltype(R), size(R)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate and return the first M moments\n\nmu_m = langle R T_m(A) R rangle\n\nif R is a vector, and if R is a matrix then\n\nmu_m = frac1N sum_n=1^N langle R_n T_m(A) R_n rangle\n\nwhere R_n rangle is the n'th column of R.\n\n\n\n\n\nkpm_moments(\n M::Int, A, bounds, U::T, V::T,\n tmp = zeros(eltype(V), size(V)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate and return the first M moments\n\nmu_m = langle U T_m(A) V rangle\n\nif U and V are vector, and if they are matrices then\n\nmu_m = frac1N sum_n=1^N langle U_n T_m(A) V_n rangle\n\nwhere U_n rangle and V_n rangle are are the n'th columns of each matrix.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_moments!","page":"API","title":"SmoQyKPMCore.kpm_moments!","text":"kpm_moments!(\n μ::AbstractVector, A, bounds, R::T,\n tmp = zeros(eltype(R), size(R)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate the moments\n\nmu_m = langle R T_m(A) R rangle\n\nif R is a vector, and if R is a matrix then\n\nmu_m = frac1N sum_n=1^N langle R_n T_m(A) R_n rangle\n\nwhere R_n rangle is the n'th column of R.\n\n\n\n\n\nkpm_moments!(\n μ::AbstractVector, A, bounds, U::T, V::T,\n tmp = zeros(eltype(V), size(V)..., 3)\n) where {T<:AbstractVecOrMat}\n\nCalculate the moments\n\nmu_m = langle U T_m(A) V rangle\n\nif U and V are vector, and if they are matrices then\n\nmu_m = frac1N sum_n=1^N langle U_n T_m(A) V_n rangle\n\nwhere U_n rangle and V_n rangle are are the n'th columns of each matrix.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_density","page":"API","title":"SmoQyKPMCore.kpm_density","text":"kpm_density(\n N::Int, μ::T, W\n) where {T<:AbstractVector}\n\nGiven the KPM moments mu, evaluate the corresponding spectral density rho(epsilon) at N energy points epsilon, where W is the spectral bandwidth or the bounds of the spectrum. This function returns both rho(epsilon) and epsilon.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_density!","page":"API","title":"SmoQyKPMCore.kpm_density!","text":"kpm_density!(\n ρ::AbstractVector{T},\n ϵ::AbstractVector{T},\n μ::AbstractVector{T},\n bounds,\n r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01)\n) where {T<:AbstractFloat}\n\nGiven the KPM moments mu, calculate the spectral density rho(epsilon) and corresponding energies epsilon where it is evaluated. Here bounds is bounds the eigenspecturm, such that the bandwidth is given by W = bounds[2] - bounds[1].\n\n\n\n\n\nkpm_density!(\n ρ::AbstractVector{T},\n ϵ::AbstractVector{T},\n μ::AbstractVector{T},\n W::T,\n r2rplan = FFTW.plan_r2r!(ρ, FFTW.REDFT01),\n) where {T<:AbstractFloat}\n\nGiven the KPM moments mu, calculate the spectral density rho(epsilon) and corresponding energies epsilon where it is evaluated. Here W is the spectral bandwidth.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_dot","page":"API","title":"SmoQyKPMCore.kpm_dot","text":"kpm_dot(\n A, coefs::AbstractVector, bounds, R::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nIf R is a single vector, then calculate the inner product\n\nbeginalign*\nS = langle R F(A) R rangle \nS = sum_m=1^M langle R c_m T_m(A^prime) R rangle\nendalign*\n\nwher A^prime is the scaled version of A. If R is a matrix, then calculate\n\nbeginalign*\nS = langle R F(A) R rangle \nS = frac1N sum_n=1^N sum_m=1^M langle R_n c_m T_m(A^prime) R_n rangle\nendalign*\n\nwhere R_n rangle is a column of R.\n\n\n\n\n\nkpm_dot(\n A, coefs::AbstractVector, bounds, U::T, V::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nIf U and V are single vectors, then calculate the inner product\n\nbeginalign*\nS = langle U F(A) V rangle \n = sum_m=1^M langle U c_m T_m(A^prime) V rangle\nendalign*\n\nwher A^prime is the scaled version of A. If U and V are matrices, then calculate\n\nbeginalign*\nS = langle U F(A) V rangle \n = frac1N sum_n=1^N sum_m=1^M langle U_n c_m T_m(A^prime) V_n rangle\nendalign*\n\nwhere U_n rangle and V_n rangle are the columns of each matrix.\n\n\n\n\n\nkpm_dot(\n A, kpm_expansion::KPMExpansion, R::T,\n tmp = zeros(eltype(R), size(R)..., 3)\n) where {T<:AbstractVecOrMat}\n\nIf R is a single vector, then calculate the inner product\n\nbeginalign*\nS = langle R F(A) R rangle \nS = sum_m=1^M langle R c_m T_m(A^prime) R rangle\nendalign*\n\nwher A^prime is the scaled version of A. If R is a matrix, then calculate\n\nbeginalign*\nS = langle R F(A) R rangle \nS = frac1N sum_n=1^N sum_m=1^M langle R_n c_m T_m(A^prime) R_n rangle\nendalign*\n\nwhere R_n rangle is a column of R.\n\n\n\n\n\nkpm_dot(\n A, kpm_expansion::KPMExpansion, U::T, V::T,\n tmp = zeros(eltype(V), size(V)..., 3)\n) where {T<:AbstractVecOrMat}\n\nIf U and V are single vectors, then calculate the inner product\n\nbeginalign*\nS = langle U F(A) V rangle \n = sum_m=1^M langle U c_m T_m(A^prime) V rangle\nendalign*\n\nwher A^prime is the scaled version of A. If U and V are matrices, then calculate\n\nbeginalign*\nS = langle U F(A) V rangle \n = frac1N sum_n=1^N sum_m=1^M langle U_n c_m T_m(A^prime) V_n rangle\nendalign*\n\nwhere U_n rangle and V_n rangle are the columns of each matrix.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_mul","page":"API","title":"SmoQyKPMCore.kpm_mul","text":"kpm_mul(\n A, coefs::AbstractVector, bounds, v::T, tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluate and return the vector v^prime = F(A) cdot v where F(A) is represented by the Chebyshev expansion. For more information refer to kpm_mul!.\n\n\n\n\n\nkpm_mul(A, kpm_expansion::KPMExpansion, v::T) where {T<:AbstractVecOrMat}\n\nEvaluate and return the vector v^prime = F(A) cdot v where F(A) is represented by the Chebyshev expansion. For more information refer to kpm_mul!.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_mul!","page":"API","title":"SmoQyKPMCore.kpm_mul!","text":"kpm_mul!(\n v′::T, A, coefs::AbstractVector, bounds, v::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluates v^prime = F(A) cdot v, writing the result to v′, where F(A) is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate u = Acdot v, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. The vector coefs contains Chebyshev expansion coefficients to approximate F(A), where the eigenspectrum of A is contained in the interval (bounds[1], bounds[2]) specified by the bounds argument. The vector v is vector getting multiplied by the Chebyshev expansion for F(A). Lastly, tmp is an array used to avoid dynamic memory allocations.\n\n\n\n\n\nkpm_mul!(\n v′::T, A, kpm_expansion::KPMExpansion, v::T,\n tmp = zeros(eltype(v), size(v)..., 3)\n) where {T<:AbstractVecOrMat}\n\nEvaluates v^prime = F(A) cdot v, writing the result to v′, where F(A) is represented by the Chebyshev expansion. Here A is either a function that can be called as A(u,v) to evaluate u = Acdot v, modifying u in-place, or is a type for which the operation mul!(u, A, v) is defined. Lastly, the array tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_eval","page":"API","title":"SmoQyKPMCore.kpm_eval","text":"kpm_eval(x::AbstractFloat, coefs, bounds)\n\nEvaluate F(x) where x is real number in the interval bounds[1] < x < bound[2], and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(A::AbstractMatrix, coefs, bounds)\n\nEvaluate and return the matrix F(A) where A is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(x::T, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}\n\nEvaluate F(x) where x is real number in the interval bounds[1] < x < bound[2], and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\nkpm_eval(A::AbstractMatrix{T}, kpm_expansion::KPMExpansion{T}) where {T<:AbstractFloat}\n\nEvaluate and return the matrix F(A) where A is an operator with strictly real eigenvalues and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.kpm_eval!","page":"API","title":"SmoQyKPMCore.kpm_eval!","text":"kpm_eval!(\n F::AbstractMatrix, A, coefs::AbstractVector, bounds,\n tmp = zeros(eltype(F), size(F)..., 3)\n)\n\nEvaluate and write the matrix F(A) to F, where A is an operator with strictly real eigenvalues that fall in the interval (bounds[1], bounds[2]) specified by the bounds argument, and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\nkpm_eval!(\n F::AbstractMatrix, A, kpm_expansion::KPMExpansion,\n tmp = zeros(eltype(F), size(F)..., 3)\n)\n\nEvaluate and write the matrix F(A) to F, where A is an operator with strictly real eigenvalues and the function F(bullet) is represented by a Chebyshev expansion with coefficients given by the vector coefs. Lastly, the array tmp is used to avoid dynamic memory allocations.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.apply_jackson_kernel","page":"API","title":"SmoQyKPMCore.apply_jackson_kernel","text":"apply_jackson_kernel(coefs)\n\nReturn the Chebyshev expansion coefficients transformed by the Jackson kernel.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.apply_jackson_kernel!","page":"API","title":"SmoQyKPMCore.apply_jackson_kernel!","text":"apply_jackson_kernel!(coefs)\n\nModify the Chebyshev expansion coefficients by applying the Jackson kernel to them.\n\n\n\n\n\napply_jackson_kernel!(kpm_expansion::KPMExpansion)\n\nModify the Chebyshev expansion coefficients by applying the Jackson kernel to them.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.lanczos","page":"API","title":"SmoQyKPMCore.lanczos","text":"lanczos(niters, v, A, S = I, rng = Random.default_rng())\n\nUse niters Lanczos iterations to find a truncated tridiagonal representation of Acdot S, up to similarity transformation. Here, A is any Hermitian matrix, while S is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, S = I. The extension to non-identity matrices S is as follows: Each matrix-vector product Acdot v becomes (A S) cdot v, and each vector inner product w^dagger cdot v becomes w^dagger cdot S cdot v. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for A and S, respectively.\n\nThis function returns a SymTridiagonal matrix. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.\n\nSimilar generalizations of Lanczos have been considered in [2] and [3].\n\n\n\n1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),\n\nhttps://doi.org/10.1093%2Fimamat%2F10.3.373.\n\n2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),\n\nhttps://doi.org/10.1090/s0025-5718-1982-0669648-0\n\n3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),\n\nhttps://doi.org/10.1016/j.commatsci.2011.02.021.\n\n\n\n\n\n","category":"function"},{"location":"api/#SmoQyKPMCore.lanczos!","page":"API","title":"SmoQyKPMCore.lanczos!","text":"lanczos!(\n αs::AbstractVector, βs::AbstractVector, v::AbstractVector,\n A, S = I,\n tmp::AbstractMatrix = zeros(eltype(v), length(v), 5);\n rng = Random.default_rng()\n)\n\nlanczos!(\n αs::AbstractMatrix, βs::AbstractMatrix, v::AbstractMatrix,\n A, S = I,\n tmp::AbstractArray = zeros(eltype(v), size(v)..., 5);\n rng = Random.default_rng()\n)\n\nUse Lanczos iterations to find a truncated tridiagonal representation of Acdot S, up to similarity transformation. Here, A is any Hermitian matrix, while S is both Hermitian and positive definite. Traditional Lanczos uses the identity matrix, S = I. The extension to non-identity matrices S is as follows: Each matrix-vector product Acdot v becomes (A S) cdot v, and each vector inner product w^dagger cdot v becomes w^dagger cdot S cdot v. The implementation below follows Wikipedia, and is the most stable of the four variants considered by Paige [1]. This implementation introduces additional vector storage so that each Lanczos iteration requires only one matrix-vector multiplication for A and S, respectively.\n\nThe number of Lanczos iterations performed equals niters = length(αs), and niters - 1 == length(βs). This function returns a SymTridiagonal matrix based on the contents of the vectors αs and βs. Note that the eigmin and eigmax routines have specialized implementations for a SymTridiagonal matrix type.\n\nNote that if αs, βs and v are all matrices, then each column of v is treated as a seperate vector, and a vector of SymTridiagonal of length size(v,2) will be returned.\n\nSimilar generalizations of Lanczos have been considered in [2] and [3].\n\n\n\n1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),\n\nhttps://doi.org/10.1093%2Fimamat%2F10.3.373.\n\n2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),\n\nhttps://doi.org/10.1090/s0025-5718-1982-0669648-0\n\n3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),\n\nhttps://doi.org/10.1016/j.commatsci.2011.02.021.\n\n\n\n\n\n","category":"function"},{"location":"usage/","page":"Usage","title":"Usage","text":"EditURL = \"../../examples/usage.jl\"","category":"page"},{"location":"usage/#Usage","page":"Usage","title":"Usage","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"Here we demonstrate the basic usage for the SmoQyKPMCore package. Let us first import the relevant packages we will want to use in this example.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"using LinearAlgebra\nusing SparseArrays","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Package for making figures","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"using CairoMakie\nCairoMakie.activate!(type = \"svg\")\n\nusing SmoQyKPMCore","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"In this usage example we will consider a 1D chain tight-binding model","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"hatH = -t sum_isigma (hatc_i+1sigma^dagger hatc_isigma + rm Hc)\n = sum_sigma hatc_isigma^dagger H_ij hatc_jsigma","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"where hatc_isigma^dagger (hatc_isigma) creates (annihilates) a spin-sigma electron on site i in the lattice, and t is the nearest-neighbor hopping amplitude.","category":"page"},{"location":"usage/#Density-Matrix-Approximation","page":"Usage","title":"Density Matrix Approximation","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"With H the Hamiltonian matrix, the corresponding density matrix is given by rho = f(H) where","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"f(epsilon) = frac11+e^beta (epsilon-mu) = frac12 left 1 + tanhleft(tfracbeta(epsilon-mu)2right) right","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"is the Fermi function, with beta = 1T the inverse temperature and mu the chemical potential. Here we will applying the kernel polynomial method (KPM) algorithm to approximate the density matrix rho by a Chebyshev polynomial expansion.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us first define the relevant model parameter values.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Nearest-neighbor hopping amplitude.\nt = 1.0\n\n# Chemical potential.\nμ = 0.0\n\n# System size.\nL = 16\n\n# Inverse temperature.\nβ = 4.0;\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now we need to construct the Hamiltonian matrix H.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"H = zeros(L,L)\nfor i in 1:L\n j = mod1(i+1,L)\n H[j,i], H[i,j] = -t, -t\nend\nH","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"We also need to define the Fermi function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Fermi function.\nfermi(ϵ, μ, β) = (1+tanh(β*(ϵ-μ)/2))/2","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us calculate the exact density matrix rho so that we may asses the accuracy of the KPM method.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Diagonalize the Hamiltonian matrix.\nϵ, U = eigen(H)\n\n# Calculate the density matrix.\nρ = U * Diagonal(fermi.(ϵ, μ, β)) * adjoint(U)","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now let use the KPMExpansion type to approximate rho. We will need to define the order M of the expansion and give approximate bounds for the eigenspectrum of H, making sure to overestimate the the true interval spanned by the eigenvalues of H. Because we are considering the simple non-interacting model, the exact eigenspectrum of H is known and spans the interval epsilon in -2t 2t","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Define eigenspectrum bounds.\nbounds = (-3.0t, 3.0t)\n\n# Define order of Chebyshev expansion used in KPM approximation.\nM = 10\n\n# Initialize expansion.\nkpm_expansion = KPMExpansion(x -> fermi(x, μ, β), bounds, M);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now let us test and see how good a job it does approximating the density matrix rho, using the kpm_eval! function to do so.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Initialize KPM density matrix approxiatmion\nρ_kpm = similar(H)\n\n# Initialize additional matrices to avoid dynamic memory allocation.\nmtmp = zeros(L,L,3)\n\n# Calculate KPM density matrix approximation.\nkpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)\n\n# Check how good an approximation it is.\nprintln(\"Matrix Error = \", norm(ρ_kpm - ρ)/norm(ρ) )","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"It is also possible to efficiently multiply vectors by KPM approximation to the density matrix rho. Let us test this functionality with the kpm_mul! on a random vector, seeing how accurate the result is.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Initialize random vector.\nv = randn(L)\n\n# Calculate exact product with density matrix.\nρv = ρ * v\n\n# Initialize vector to contain approximate product with densith matrix.\nρv_kpm = similar(v)\n\n# Initialize to avoid dynamic memory allocation.\nvtmp = zeros(L, 3)\n\n# Calculate approximate product with density matrix.\nkpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)\n\n# Check how good the approximation is.\nprintln(\"Vector Error = \", norm(ρv_kpm - ρv) / norm(ρv) )","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us now provide the KPM approximation with better bounds on the eigenspectrum and also increase the order the of the expansion and see how the result improves. We will use the kpm_update! function to update the KPMExpansion in-place.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Define eigenspectrum bounds.\nbounds = (-2.5t, 2.5t)\n\n# Define order of Chebyshev expansion used in KPM approximation.\nM = 100\n\n# Update KPM approximation.\nkpm_update!(kpm_expansion, x -> fermi(x, μ, β), bounds, M)\n\n# Calculate KPM density matrix approximation.\nkpm_eval!(ρ_kpm, H, kpm_expansion, mtmp)\n\n# Check how good an approximation it is.\nprintln(\"Matrix Error = \", norm(ρ_kpm - ρ)/norm(ρ) )\n\n# Calculate approximate product with density matrix.\nkpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp)\n\n# Check how good the approximation is.\nprintln(\"Vector Error = \", norm(ρv_kpm - ρv) / norm(ρv) )","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Now let me quickly demonstrate how we may approximate the the trace rho using a set of random vector R_n using the kpm_dot function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Calculate exact trace.\ntrρ = tr(ρ)\nprintln(\"Exact trace = \", trρ)\n\n# Number of random vectors\nN = 100\n\n# Initialize random vectors.\nR = randn(L, N)\n\n# Initialize array to avoid dynamic memory allocation\nRtmp = zeros(L, N, 3)\n\n# Approximate trace of density matrix.\ntrρ_approx = kpm_dot(H, kpm_expansion, R, Rtmp)\nprintln(\"Approximate trace = \", trρ_approx)\n\n# Report the error in the approximation.\nprintln(\"Trace esitimate error = \", abs(trρ_approx - trρ)/trρ)","category":"page"},{"location":"usage/#Density-of-States-Approximation","page":"Usage","title":"Density of States Approximation","text":"","category":"section"},{"location":"usage/","page":"Usage","title":"Usage","text":"Next let us demonstrate how we may approximate the density of states mathcalN(epsilon) for a 1D chain tight-binding model. The first step a very larger Hamiltonian matrix H which we represent as a spare matrix.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Size of 1D chain considered\nL = 10_000\n\n# Construct sparse Hamiltonian matrix.\nrows = Int[]\ncols = Int[]\nfor i in 1:L\n j = mod1(i+1, L)\n push!(rows,i)\n push!(cols,j)\n push!(rows,j)\n push!(cols,i)\nend\nvals = fill(-t, length(rows))\nH = sparse(rows, cols, vals)","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Let us now calculate the fist M moments mu_m using N random vectors using the function kpm_moments.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Number of moments to calculate.\nM = 250\n\n# Number of random vectors used to approximate moments.\nN = 100\n\n# Initialize random vectors.\nR = randn(L, N)\n\n# Calculate the moments.\nμ = kpm_moments(M, H, bounds, R);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Having calculate the moments, let us next evaluate the density of states mathcalN(epsilon) at P points with the kpm_density function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Number of points at which to evaluate the density of states.\nP = 1000\n\n# Evaluate density of states.\ndos, ϵ = kpm_density(P, μ, bounds);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Without regularization, the approximate for the density of states generated above will have clearly visible Gibbs oscillations. To partially suppress these artifacts, we apply the Jackson kernel to the moments mu using the apply_jackson_kernel function.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"# Apply Jackson kernel.\nμ_jackson = apply_jackson_kernel(μ)\n\n# Evaluate density of states.\ndos_jackson, ϵ_jackson = kpm_density(P, μ_jackson, bounds);\nnothing #hide","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"Have approximated the density of states, let us now plot the result.","category":"page"},{"location":"usage/","page":"Usage","title":"Usage","text":"fig = Figure(\n size = (700, 400),\n fonts = (; regular= \"CMU Serif\"),\n figure_padding = 10\n)\n\nax = Axis(\n fig[1, 1],\n aspect = 7/4,\n xlabel = L\"\\epsilon\", ylabel = L\"\\mathcal{N}(\\epsilon)\",\n xticks = range(start = bounds[1], stop = bounds[2], length = 5),\n xlabelsize = 30, ylabelsize = 30,\n xticklabelsize = 24, yticklabelsize = 24,\n)\n\nlines!(\n ϵ, dos,\n linewidth = 2.0,\n alpha = 1.0,\n color = :red,\n linestyle = :solid,\n label = \"Dirichlet\"\n)\n\nlines!(\n ϵ_jackson, dos_jackson,\n linewidth = 3.0,\n alpha = 1.0,\n color = :black,\n linestyle = :solid,\n label = \"Jackson\"\n)\n\nxlims!(\n ax, bounds[1], bounds[2]\n)\n\nylims!(\n ax, 0.0, 1.05 * maximum(dos)\n)\n\naxislegend(\n ax, halign = :center, valign = :top, labelsize = 30\n)\n\nfig","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = SmoQyKPMCore","category":"page"},{"location":"#SmoQyKPMCore","page":"Home","title":"SmoQyKPMCore","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Documentation for SmoQyKPMCore package. The SmoQyKPMCore package implements and exports an optimized, low-level implementation of the Kernel Polynomial Method (KPM) algorithm for approximating functions of operators with strictly real, bounded eigenvalues via a Chebyshev polynomial expansion.","category":"page"},{"location":"#Funding","page":"Home","title":"Funding","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"The development of this code was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0022311.","category":"page"}] } diff --git a/dev/usage/86e74be1.svg b/dev/usage/bb1f32c8.svg similarity index 54% rename from dev/usage/86e74be1.svg rename to dev/usage/bb1f32c8.svg index dd35e62..75d396e 100644 --- a/dev/usage/86e74be1.svg +++ b/dev/usage/bb1f32c8.svg @@ -2,278 +2,278 @@ - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + - + - - - - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - - + + - + - - + + - + - + - - + + - - + + - + - + - + - + - + - + - + - + - + - + - + - + - + - - + + - + - + - + - + - + - + - + - - + + - + - + - + - + - + - + - + - + - + - + - - + + - + - + - + - + - + - + - + - + - + - + - + - + @@ -285,261 +285,261 @@ - - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + @@ -548,59 +548,59 @@ - - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + diff --git a/dev/usage/index.html b/dev/usage/index.html index 6151663..9426ed4 100644 --- a/dev/usage/index.html +++ b/dev/usage/index.html @@ -89,7 +89,7 @@ kpm_mul!(ρv_kpm, H, kpm_expansion, v, vtmp) # Check how good the approximation is. -println("Vector Error = ", norm(ρv_kpm - ρv) / norm(ρv) )
Vector Error = 0.015238739082477082

Let us now provide the KPM approximation with better bounds on the eigenspectrum and also increase the order the of the expansion and see how the result improves. We will use the kpm_update! function to update the KPMExpansion in-place.

# Define eigenspectrum bounds.
+println("Vector Error = ", norm(ρv_kpm - ρv) / norm(ρv) )
Vector Error = 0.013942676173184827

Let us now provide the KPM approximation with better bounds on the eigenspectrum and also increase the order the of the expansion and see how the result improves. We will use the kpm_update! function to update the KPMExpansion in-place.

# Define eigenspectrum bounds.
 bounds = (-2.5t, 2.5t)
 
 # Define order of Chebyshev expansion used in KPM approximation.
@@ -109,7 +109,27 @@
 
 # Check how good the approximation is.
 println("Vector Error = ", norm(ρv_kpm - ρv) / norm(ρv) )
Matrix Error = 1.031960839804449e-14
-Vector Error = 1.127015333309712e-14

Density of States Approximation

Next let us demonstrate how we may approximate the density of states $\mathcal{N}(\epsilon)$ for a 1D chain tight-binding model. The first step a very larger Hamiltonian matrix $H,$ which we represent as a spare matrix.

# Size of 1D chain considered
+Vector Error = 9.399834397240509e-15

Now let me quickly demonstrate how we may approximate the the trace $\rho$ using a set of random vector $R_n$ using the kpm_dot function.

# Calculate exact trace.
+trρ = tr(ρ)
+println("Exact trace = ", trρ)
+
+# Number of random vectors
+N = 100
+
+# Initialize random vectors.
+R = randn(L, N)
+
+# Initialize array to avoid dynamic memory allocation
+Rtmp = zeros(L, N, 3)
+
+# Approximate trace of density matrix.
+trρ_approx = kpm_dot(H, kpm_expansion, R, Rtmp)
+println("Approximate trace = ", trρ_approx)
+
+# Report the error in the approximation.
+println("Trace esitimate error = ", abs(trρ_approx - trρ)/trρ)
Exact trace = 8.000000000000002
+Approximate trace = 7.91420719844715
+Trace esitimate error = 0.010724100194106521

Density of States Approximation

Next let us demonstrate how we may approximate the density of states $\mathcal{N}(\epsilon)$ for a 1D chain tight-binding model. The first step a very larger Hamiltonian matrix $H,$ which we represent as a spare matrix.

# Size of 1D chain considered
 L = 10_000
 
 # Construct sparse Hamiltonian matrix.
@@ -206,4 +226,4 @@
     ax, halign = :center, valign = :top, labelsize = 30
 )
 
-fig
Example block output +figExample block output