This function allows the user to build an array of matrices, where each matrix is a symmetric matrix of correlators at a given timeslice.
It takes as input:
- `diag::Vector{Array}`: vector of correlators. Each correlator will constitute a diagonal element of the matrices A[i] where i runs over the timeslices, i.e.
`A[i][1,1] = diag[1], .... A[i][n,n] = diag[n]`
`Given n=length(diag)`, the matrices will have dimension n*n
- `upper::Vector{Array}`: vector of correlators liying on the upper diagonal position.
This method computes the energy level from the eigenvalues according to:
``E_i(t) = \log(λ(t) / λ(t+1))``
where `i=1,..,n` with `n=length(evals[1])` and `t=1,..,T` total time slices.
It returns a vector array en where each entry en[i][t] contains the i-th states energy at time t
Examples:
```@example
## load data
pp_data = read_mesons(path, "G5", "G5")
pa_data = read_mesons(path, "G5", "G0G5")
aa_data = read_mesons(path, "G0G5", "G0G5")
## create Corr struct
corr_pp = corr_obs.(pp_data)
corr_pa = corr_obs.(pa_data)
corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
## set up matrices
corr_diag = [corr_pp[1], corr_aa[1]]
corr_upper = [corr_pa[1]]
matrices = [corr_diag, corr_upper]
## solve the GEVP
evals = getall_eigvals(matrices, 5) #where t_0=5
en = energies(evals)
Julia> en[i] # i-th state energy at each timeslice
```
"""
function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}}};wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}},Nothing}=nothing)
time=length(evals)
n=length(evals[1])
eff_en=Array{Array{uwreal}}(undef,n)
aux_en=Vector{uwreal}(undef,time-1)
foriin1:n
fortin1:time-1
ratio=evals[t][i]/evals[t+1][i]
aux_en[t]=0.5*log(ratio*ratio)
end
#uwerr.(aux_en)
#isnothing(wpm) ? uwerr.(aux_en) : [uwerr(a, wpm) for a in aux_en]
eff_en[i]=copy(aux_en)
end
returneff_en
end
@docraw"""
getall_eigvals(a::Vector{Matrix}, t0; iter=30 )
This function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one
at index t0, i.e:
``C(t_i)v_i = λ_i C(t_0) v_i``, with i=1:lenght(a)
It takes as input:
- `a::Vector{Matrix}` : a vector of matrices
- `t0::Int64` : idex value at which the fixed matrix is taken
- `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues
It returns:
- `res` = Vector{Vector{uwreal}}
where `res[i]` are the generalised eigenvalues of the i-th matrix of the input array.
Examples:
```@example
## load data
pp_data = read_mesons(path, "G5", "G5")
pa_data = read_mesons(path, "G5", "G0G5")
aa_data = read_mesons(path, "G0G5", "G0G5")
## create Corr struct
corr_pp = corr_obs.(pp_data)
corr_pa = corr_obs.(pa_data)
corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
## set up matrices
corr_diag = [corr_pp[1], corr_aa[1]]
corr_upper = [corr_pa[1]]
matrices = [corr_diag, corr_upper]
## solve the GEVP
#evals = getall_eigvals(matrices, 5) #where t_0=5
Julia>
```
"""
function getall_eigvals(a::Vector{Matrix{uwreal}},t0::Int64;iter=5)
n=length(a)
res=Vector{Vector{uwreal}}(undef,n)
[res[i]=uweigvals(a[i],a[t0],iter=iter)fori=1:n]
returnres
end
@docraw"""
uweigvals(a::Matrix{uwreal}; iter = 30)
uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
This function computes the eigenvalues of a matrix of uwreal objects.
If a second matrix b is given as input, it returns the generalised eigenvalues instead.
It takes as input:
- `a::Matrix{uwreal}` : a matrix of uwreal
- `b::Matrix{uwreal}` : a matrix of uwreal, optional
- `iter=30`: optional flag to set the iterations of the qr algorithm used to solve the eigenvalue problem
It returns:
- `res = Vector{uwreal}`: a vector where each elements is an eigenvalue
```@example
a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
res = uweigvals(a) ##eigenvalues
res1 = uweigvals(a,b) ## generalised eigenvalues
```
"""
function uweigvals(a::Matrix{uwreal};iter::Int64=5,diag::Bool=true)
This method uses backward propagation to compute the inverse of an upper triangular matrix u. Note that the inverse of u must be upper triangular as well.
Tests executed so far always confirmed this property.
This function allows the user to build an array of matrices, where each matrix is a symmetric matrix of correlators at a given timeslice.
It takes as input:
- `diag::Vector{Array}`: vector of correlators. Each correlator will constitute a diagonal element of the matrices A[i] where i runs over the timeslices, i.e.
`A[i][1,1] = diag[1], .... A[i][n,n] = diag[n]`
`Given n=length(diag)`, the matrices will have dimension n*n
- `upper::Vector{Array}`: vector of correlators liying on the upper diagonal position.
This method computes the energy level from the eigenvalues according to:
``E_i(t) = \log(λ(t) / λ(t+1))``
where `i=1,..,n` with `n=length(evals[1])` and `t=1,..,T` total time slices.
It returns a vector array en where each entry en[i][t] contains the i-th states energy at time t
Examples:
```@example
## load data
pp_data = read_mesons(path, "G5", "G5")
pa_data = read_mesons(path, "G5", "G0G5")
aa_data = read_mesons(path, "G0G5", "G0G5")
## create Corr struct
corr_pp = corr_obs.(pp_data)
corr_pa = corr_obs.(pa_data)
corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
## set up matrices
corr_diag = [corr_pp[1], corr_aa[1]]
corr_upper = [corr_pa[1]]
matrices = [corr_diag, corr_upper]
## solve the GEVP
evals = getall_eigvals(matrices, 5) #where t_0=5
en = energies(evals)
Julia> en[i] # i-th state energy at each timeslice
```
"""
function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}}};wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}},Nothing}=nothing)
time=length(evals)
n=length(evals[1])
eff_en=Array{Array{uwreal}}(undef,n)
aux_en=Vector{uwreal}(undef,time-1)
foriin1:n
fortin1:time-1
ratio=evals[t][i]/evals[t+1][i]
aux_en[t]=0.5*log(ratio*ratio)
end
#uwerr.(aux_en)
#isnothing(wpm) ? uwerr.(aux_en) : [uwerr(a, wpm) for a in aux_en]
eff_en[i]=copy(aux_en)
end
returneff_en
end
@docraw"""
getall_eigvals(a::Vector{Matrix}, t0; iter=30 )
This function solves a GEVP problem, returning the eigenvalues, for a list of matrices, taking as generalised matrix the one
at index t0, i.e:
``C(t_i)v_i = λ_i C(t_0) v_i``, with i=1:lenght(a)
It takes as input:
- `a::Vector{Matrix}` : a vector of matrices
- `t0::Int64` : idex value at which the fixed matrix is taken
- `iter=30` : the number of iterations of the qr algorithm used to extract the eigenvalues
It returns:
- `res` = Vector{Vector{uwreal}}
where `res[i]` are the generalised eigenvalues of the i-th matrix of the input array.
Examples:
```@example
## load data
pp_data = read_mesons(path, "G5", "G5")
pa_data = read_mesons(path, "G5", "G0G5")
aa_data = read_mesons(path, "G0G5", "G0G5")
## create Corr struct
corr_pp = corr_obs.(pp_data)
corr_pa = corr_obs.(pa_data)
corr_aa = corr_obs.(aa_data) # array of correlators for different \mu_q combinations
## set up matrices
corr_diag = [corr_pp[1], corr_aa[1]]
corr_upper = [corr_pa[1]]
matrices = [corr_diag, corr_upper]
## solve the GEVP
#evals = getall_eigvals(matrices, 5) #where t_0=5
Julia>
```
"""
function getall_eigvals(a::Vector{Matrix{uwreal}},t0::Int64;iter=5)
n=length(a)
res=Vector{Vector{uwreal}}(undef,n)
[res[i]=uweigvals(a[i],a[t0],iter=iter)fori=1:n]
returnres
end
@docraw"""
uweigvals(a::Matrix{uwreal}; iter = 30)
uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
This function computes the eigenvalues of a matrix of uwreal objects.
If a second matrix b is given as input, it returns the generalised eigenvalues instead.
It takes as input:
- `a::Matrix{uwreal}` : a matrix of uwreal
- `b::Matrix{uwreal}` : a matrix of uwreal, optional
- `iter=30`: optional flag to set the iterations of the qr algorithm used to solve the eigenvalue problem
It returns:
- `res = Vector{uwreal}`: a vector where each elements is an eigenvalue
```@example
a = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
b = Matrix{uwreal}(nothing, n,n) ## n*n matrix of uwreal with nothing entries
res = uweigvals(a) ##eigenvalues
res1 = uweigvals(a,b) ## generalised eigenvalues
```
"""
function uweigvals(a::Matrix{uwreal};iter::Int64=5,diag::Bool=true)
This method uses backward propagation to compute the inverse of an upper triangular matrix u. Note that the inverse of u must be upper triangular as well.
Tests executed so far always confirmed this property.
