Commit bb15ee61 by Alessandro

Introducing new functions and documentations for eigenvalue and eigenvector routines

parent 62815992
 #=using ForwardDiff, LinearAlgebra #= using ForwardDiff, LinearAlgebra for op in (:eigvals, :eigvecs) @eval function LinearAlgebra.\$op(a::Matrix{uwreal}) ... ... @@ -51,7 +52,6 @@ Julia> 4*4 Array{uwreal,2} ``` """ function get_matrix(corr_diag::Vector{Vector{uwreal}}, corr_upper::Vector{Vector{uwreal}}) println("a") time = length(corr_diag[1]) # total time slices n = length(corr_diag) # matrix dimension d = length(corr_upper) # upper elements ... ... @@ -103,81 +103,46 @@ function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }) end @doc raw""" uwgevp_tot(mat_list::Vector{Array}, tnot::Int64; iter = 30, evec = false) Given a list of uwreal matrices and an integer value tnot, this method solves the generalised eigenvalue problem for all the matrices in the list, using the same auxiliary matrix mat_list[tnot]. By default, it returns a vector evals where each entry evals[i] is a vector of uwreal containing the generalised eigenvalues of the matrix mat_list[i]. If evec is set to true, the methods returns also an array evecs where each entry evecs[i] is a matrix of uwreal. The columns of such matrix are the eigenvectors associated with eigenvalues evals[i]. The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping. """ function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, evec::Bool = false) n = length(mat_list) evals = Array{Array{uwreal}}(undef, 0 ) if !evec c_inv = invert(mat_list[tnot]) for i =1:n aux_evals = uwgevp(mat_list[i], c_inv, iter = iter) push!(evals , get_diag(aux_evals)) end return evals else evecs = Array{Matrix}(undef, 0) for i = 2:n-4 c_inv = invert(mat_list[i]) aux_evals, aux_evec = uwgevp(mat_list[i+3], c_inv, evec=evec, iter=iter) push!(evals, get_diag(aux_evals)) ptilda = norm_evec(aux_evec, mat_list[i]) m = uwdot(mat_list[i], ptilda) push!(evecs, m) end return evals, evecs end end function uwgevp_tot(mat_list::Vector{Matrix}; iter::Int64=30, delta_t::Int64=3) n = length(mat_list) - delta_t evecs = Array{Matrix}(undef, 0) for i = 2:n c_inv = invert(mat_list[i]) usless, temp_evec = uwgevp(mat_list[i+delta_t], c_inv, evec=true, iter=iter) push!(evecs, uwdot(mat_list[i], norm_evec(temp_evec, mat_list[i]))) end return evecs end """ uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter = 30, evec = false) Given two matrices of uwreal C(t) and C(t0)^-1, this method solves the generalized eigenvalue problem C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matrix of eigenvalues. If evec = true also the eigenvectors are returned as columns of the orthogonal matrix u. The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping. """ function uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::Matrix{uwreal}; iter::Int64 = 30, evec::Bool = false) c = uwdot(c_tnot_inv, c_t) #matrix to diagonalize return uwevp(c, iter = iter, evec = evec) end """ uwevp(a::Matrix{uwreal}; iter = 30, evec = false) Return the eigenvalues of the uwreal matrix a. If evec = true also the eigenvectors are returned as columns of the orthogonal matrix u. The keyword iter, set by default to 30, selects the number of iterations of the qr algorithm before stopping. """ function getall_eigvals(a::Vector{Matrix}, t0; iter=30 ) 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 mat_array = get_matrix(diag, upper_diag) evals = getall_eigvals(mat_array, 5) ''' """ function getall_eigvals(a::Vector{Matrix}, t0::Int64; iter=30 ) n = length(a) res = Vector{Vector{uwreal}}(undef, n) [res[i] = uweigvals(a[i], a[t0]) for i=1:n] return res end @doc raw""" 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 It returns: - res = Vector{uwreal}: a vector where each elements is an eigenvalue """ function uweigvals(a::Matrix{uwreal}; iter = 30) n = size(a,1) for k in 1:iter ... ... @@ -195,6 +160,26 @@ function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30) return get_diag(c) end @doc raw""" getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 ) This function solves a GEVP problem, returning the eigenvectors, for a list of matrices. C(t_i)v_i = λ_i C(t_i-delta_t) v_i, with i=1:lenght(a) Here delta_t is the time shift within the two matrices of the problem, and is kept fixed. It takes as input: - a::Vector{Matrix} : a vector of matrices - delta_t::Int64 : the fixed time shift t-t_0 - iter=30 : the number of iterations of the qr algorithm used to extract the eigenvalues It returns: - res = Vector{Matrix{uwreal}} where each res[i] is a matrix with the eigenvectors as columns Examples: '''@example mat_array = get_matrix(diag, upper_diag) evecs = getall_eigvecs(mat_array, 5) ''' """ function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30) n = length(a)-delta_t res = Vector{Matrix{uwreal}}(undef, n) ... ... @@ -202,6 +187,20 @@ function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30) return res end @doc raw""" uweigvecs(a::Matrix{uwreal}; iter = 30) uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30) This function computes the eigenvectors of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvectors instead. It takes as input: - a::Matrix{uwreal} : a matrix of uwreal - b::Matrix{uwreal} : a matrix of uwreal, optional It returns: - res = Matrix{uwreal}: a matrix where each column is an eigenvector """ function uweigvecs(a::Matrix{uwreal}; iter = 30) n = size(a,1) u = idty(n) ... ... @@ -224,6 +223,19 @@ function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30) return u end @doc raw""" uweigen(a::Matrix{uwreal}; iter = 30) uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30) This function computes the eigenvalues and the eigenvectors of a matrix of uwreal objects. If a second matrix b is given as input, it returns the generalised eigenvalues and eigenvectors instead. It takes as input: - a::Matrix{uwreal} : a matrix of uwreal - b::Matrix{uwreal} : a matrix of uwreal, optional It returns: - evals = Vector{uwreal}: a vector where each elements is an eigenvalue -evecs = Matrix{uwreal}: a matrix where the i-th column is the eigenvector of the i-th eigenvalue """ function uweigen(a::Matrix{uwreal}; iter = 30) return uweigvals(a, iter=iter), uweigvecs(a, iter=iter) end ... ... @@ -238,9 +250,6 @@ function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal}) aux_norm = (uwdot(evec[:,i], uwdot(ctnot, evec[:,i])).^2).^0.25 res_evec[:,i] = 1 ./aux_norm .* evec[:,i] end #println("res_evec is ", res_evec) #println("evec is ", evec) #println(res_evec) return res_evec end function get_diag(a::Matrix{uwreal}) ... ...
 ... ... @@ -72,7 +72,7 @@ end @doc raw""" read_mesons(path::String, g1::Union{String, Nothing}=nothing, g2::Union{String, Nothing}=nothing; id::Union{Int64, Nothing}=nothing) This faction read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures. This function read a mesons dat file at a given path and returns a vector of CData structures for different masses and Dirac structures. Dirac structures g1 and/or g2 can be passed as string arguments in order to filter correaltors. ADerrors id can be specified as argument. If is not specified, the id is fixed according to the ensemble name (example: "H400"-> id = 400) ... ...
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!