Skip to content

J
juobs
Commit 6a9891c7 authored by Alessandro 's avatar Alessandro
Browse files
Options

uwLinAlg deleted and performance improvement in juobs_linalg.jl

parent dbb425c4
Hide whitespace changes
Inline Side-by-side
Showing with 4 additions and 481 deletions
src/juobs_linalg.jl
View file @ 6a9891c7
......@@ -117,12 +117,13 @@ The columns of such matrix are the eigenvectors associated with eigenvalues eval
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)
c_inv = invert(mat_list[tnot])
if !evec
aux_evals = uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter)
aux_evals = uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter)
evals = get_diag.(aux_evals)
return evals
else
aux_evals, aux_evec = uwgevp.(mat_list, c_tnot=mat_list[tnot], iter = iter, evec= true)
aux_evals, aux_evec = uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter, evec= true)
evals = get_diag.(aux_evals)
evecs = aux_evec
return evals, evecs
......@@ -139,8 +140,7 @@ C(t)V = C(t_0)D V, where V is the eigenvector matrix and D is the diagonal matri
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::Matrix{uwreal}, iter::Int64 = 30, evec::Bool = false)
c_tnot_inv = invert(c_tnot)
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
......
src/uwLinAlg.jl deleted 100644 → 0
View file @ dbb425c4
##################################################################
# Alessandro Conigli wrote this script.
# It implements multiple linear algebra utilites optimised to work
# with uwreal data types.
# The major features so fare implemented are:
# - matrix matrix and vector vector multiplication
# - norm of uwreal vectors
# - copysign function for uwreal variables
# - creation of identity matrix of uwreal
# - element wise multiplication between uwreal variables and matrix of uwreal
#
#
##################################################################
## Setting up packages
using PyPlot, ADerrors, LsqFit, Distributions, Statistics, DelimitedFiles, ForwardDiff, LinearAlgebra
"""
get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
This method returns an array of dim T where each element is a symmetrix matrix of dimension n of uwreal correlators
at fixed time i=1..T. It takes as input:
corr_diag : vector of dimension n of correlators liying on the diagonal
corr_upper : vector of correlators liying on the upper diagonal.
Each correlator is an vector of uwreal variables of dimension T.
Example:
for i in 1:n
a[i,i] = Vector{uwreal} # vector of uwreal variables of dimension T. They will constitute the diagonal elements of the matrices
for i in 1:n-1
for j in i+1:n
a[i,j] = Vector{uwreal} # vector of uwreal variables of dimension T. They will constitute the upper diagonal elements of the matrices. A matrix
of dimension n*n has n(n-1)/2 upper diagonal elements.
Assume n=4
diagonal = Vector{Array}()
push!(diagonal, a[1,1],a[2,2],a[3,3],a[4,4])
upsize = Vector{Array}()
push!(upsize, a[1,2], a[1,3], a[1,4], a[2,3], a[2,4], a[3,4])
array_of_matrices = get_matrix(diagonal, upsize)
Julia> T-element Array{Array,1}
size(array_of_matrices)
Julia> (T,)
array_of_matrices[t] # t in 1:T
Julia> 4*4 Array{uwreal,2}
"""
function get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
time = length(corr_diag[1]) # total time slices
n = length(corr_diag) # matrix dimension
d = length(corr_upper) # upper elements
res = Vector{Array}() # array with all the matrices at each time step.
for t in 1:time
aux = Matrix{uwreal}(undef, n, n)
count=0
testing = corr_upper[1][1]
for i in range(n-1,1, step=-1)
for j in range(n, i+1, step=-1)
aux[i,j] = corr_upper[d-count][t]
aux[j,i] = corr_upper[d-count][t]
count +=1
end
aux[i,i] = corr_diag[i][t]
end
aux[n,n] = corr_diag[n][t]
push!(res, aux)
end
return res
end
"""
energies(evals::Vector{Array})
Given a vector where each entry evals[t] is a uwreal array of eigenvalues, this method computes the effective energies of the first N states, where N=dim(evals[t]).
The index t here runs from 1:T=lenght(evals), while the index i stands for the number of energy levels computed: i = length(evals[t])
It returns a vector array eff_en where each entry eff_en[t] contains the first N states energies as uwreal objects
"""
function energies(evals::Vector{Array})
time = length(evals)
n = length(evals[1])
eff_en = Vector{Array}()
for i in 1:n
aux_en = Vector{uwreal}()
for t in 1:time-1
ratio = evals[t][i] / evals[t+1][i]
push!(aux_en, log(ratio*uwcopysign(uwreal(1.0),ratio)))
end
uwerr.(aux_en)
push!(eff_en, aux_en)
end
return eff_en
end
"""
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{Array}, tnot::Int64; iter = 30, evec = false)
n = length(mat_list)
evals = Vector{Array}()
if evec == false
for i in 1:n
aux_evals = uwgevp(mat_list[i], mat_list[tnot], iter = iter)
push!(evals, get_diag(aux_evals))
end
return evals
else
evecs = Vector{Array}()
for i in 1:n
aux_evals, aux_evec= uwgevp(mat_list[i], mat_list[tnot], iter = iter, evec = true)
push!(evals, get_diag(aux_evals))
push!(evecs, aux_evec)
end
return evals, evecs
end
end
"""
uwgevp(c_t:: Matrix{uwreal}, c_tnot::Matrix{uwreal}; iter = 30, evec = false)
Given two matrices of uwreal C(t) and C(t0), 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::Matrix{uwreal}; iter = 30, evec = false)
c_tnot_inv = invert(c_tnot)
c = uwdot(c_tnot_inv, c_t) #matrix to diagonalize
if evec == false
return uwevp(c, iter = iter)
else
return uwevp(c, iter = iter, evec = true)
end
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 uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
n = size(a,1)
if evec == true
u = idty(n)
for k in 1:iter
q_k, r_k = qr(a)
a = uwdot(r_k, q_k)
u = uwdot(u, q_k)
end
return a, u
else
for k in 1:iter
q_k, r_k = qr(a)
a = uwdot(r_k, q_k)
end
return a
end
end
function get_diag(a::Matrix{uwreal})
n = size(a,1)
res = Vector{uwreal}()
for i in 1:n
push!(res, a[i,i])
end
return res
end
"""
Given an input matrix of uwreals, this function returns the inverse of that matrix using QR decomposition and
backward substitution.
"""
function invert(a::Matrix{uwreal})
#q,r = qr(hess_red(a))
q,r = qr(a)
r_inv = upper_inversion(r)
return uwdot(r_inv, transpose(q))
end
"""
qr(A::Matrix{uwreal})
qr(A::Matrix{Float64})
This methods returns the QR decomposition of a matrix A of either uwreal or Float64 variables
"""
function qr(A::Matrix{uwreal})
m,n = size(A)
Q = idty(m)
for i in 1:(n-(m==n))
H = idty(m)
H[i:m,i:m] = make_householder(A[i:m,i])
Q = uwdot(Q,H)
A = uwdot(H,A)
end
for i in 2:n
for j in 1:i-1
A[i,j] = 0.0
end
end
return Q,A
end
function qr(A::Matrix{Float64})
m,n = size(A)
Q = Matrix{Float64}(I,m,m)
for i in 1:(n-(m==n))
H = Matrix{Float64}(I,m,m)
H[i:m,i:m] = make_householder(A[i:m,i])
X = similar(Q)
Y = similar(A)
Q = mul!(X,Q,H)
A = mul!(Y,H,A)
end
return Q,A
end
"""
make_householder(a::Vector{uwreal})
make_householder(a::Vector{Float64})
This method returns the householder matrix used in qr decomposition to get an upper triangular matrix
It takes as input either a vector of uwreal or a vector of Float64
"""
function make_householder(a::Vector{uwreal})
n = length(a)
v = Vector{uwreal}(undef, n)
for i in 1:n
v[i] = a[i] / (a[1] + uwcopysign(uwnorm(a), a[1]))
end
v[1] = uwreal(1.0)
H = idty(n)
H = H - (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), reshape(v, 1 ,:))
#H = H - (2.0 / dot(value.(v),value.(v))) * uwdot(reshape(v, :,1), reshape(v, 1 ,:))
for i in 1:n
for j in 1:n
#println(H[i,j].ids)
if H[i,j]!= 0 && H[i,j]!=1
#uwerr(H[i,j])
end
end
end
return H
end
function make_householder(a::Vector{Float64})
n = length(a)
v = a / (a[1] + copysign(norm(a), a[1]))
v[1] = 1.0
H = Matrix{Float64}(I, n, n)
println(typeof(v))
y = similar(H)
println(typeof(y))
H = H .- (2.0 / dot(v,v)) * mul!( y, reshape(v, :,1), reshape(v, 1 ,:))
return H
end
"""
hess_red(a::Matrix{uwreal})
Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix.
"""
function hess_red(a::Matrix{uwreal})
n = size(a, 1)
res = Matrix{uwreal}(undef, n, n)
res = a
for j in 1:n-2
mi = idty(n)
mi_inv = idty(n)
for k in j+2:n
mi[k,j+1] = res[k,j] / res[j+1, j]
mi_inv[k,j+1] = -mi[k,j+1]
end
res = uwdot(mi_inv, uwdot(res, mi) )
end
return res
end
"""
upper_inversion(u::Matrix{uwreal})
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.
"""
function upper_inversion(u::Matrix{uwreal})
n = size(u,1)
u_inv = Matrix{uwreal}(undef, n, n)
for i in 1:n
e_i = canonical_basis(n, i)
x = backward_sub(u, e_i)
u_inv[:,i] = x
end
return u_inv
end
"""
backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
This method operates backward substitution to solve a liner system Ux = y where U is an upper triangular uwreal matrix
and y is a uwreal vector. It returns a vector of uwreal x that might be used to inverte the upper triangular matrix U.
"""
function backward_sub(u::Matrix{uwreal}, y::Vector{uwreal})
n = length(y)
x = Vector{uwreal}(undef, n)
x[n] = y[n] / u[n,n]
for i in range(n-1,1, step=-1)
temp = 0.0
for j in i+1:n
temp += u[i,j] * x[j]
end
x[i] = (y[i] - temp) /u[i,i]
if x[i] !=0
#uwerr(x[i])
end
end
return x
end
"""
canonical_basis(n::Int64, k::Int64)
Given the dimension n and the position k, canonical_basis returns a vector of uwreal with all entries set to zero
except for the k-th entry which is set to 1.0.
This method is used as input vector in backward_sub to compute the inverse of an upper triangular matrix.
"""
function canonical_basis(n::Int64, k::Int64)
e = Vector{uwreal}(undef,n)
for i in 1:n
e[i] = 0.0
end
e[k] = 1.0
return e
end
"""
uwdot(a::Vector{uwreal}, b::Vector{uwreal})
uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
Depending on the arguments, uwdot returns the vector vector or the matrix matrix product
of uwreal data type.
"""
function uwdot(a::Vector{uwreal}, b::Vector{uwreal})
n = length(a)
c = uwreal(0) #Vector{uwreal}(undef, n)
for i in 1:n
c += a[i] * b[i]
end
#uwerr(c)
return c
end
function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
n,m = size(a)
k,p = size(b)
c = Matrix{uwreal}(undef, n,p)
for i in 1:n
for j in 1:p
temp = 0.0 # uwreal()
for l in 1:m
temp += a[i,l] * b[l,j]
end
c[i,j] = temp
if c[i,j] !=0 && c[i,j]!=1
#uwerr(c[i,j])
end
end
end
return c
end
"""
uwnorm(a::Vector{uwreal})
It returns the norm of a vector of uwreal
"""
function uwnorm(a::Vector{uwreal})
norm = sqrt(uwdot(a,a))
uwerr(norm)
return norm
end
"""
Return the transpose of a matrix of uwreals
"""
function Base.transpose(a::Matrix{uwreal})
res = similar(a)
n = size(a, 1 )
for i in 1:n-1
for j in i+1:n
temp = a[i,j]
res[i,j] = a[j,i]
res[j,i] = temp
end
res[i,i] = a[i,i]
end
res[n,n] = a[n,n]
return res
end
"""
uwcopysign(a::uwreal, b::uwreal)
It implements the copysign function for uwreal variables
"""
function uwcopysign(a::uwreal, b::uwreal)
aux = copysign(value(a), value(b))
a.mean = aux
return a
end
"""
idty(n::Int64)
It returns an indetity matrix of uwreal variables given the size n
"""
function idty(n::Int64)
res = Matrix{uwreal}(undef,n,n)
for i in 1:n
for j in 1:n
if i==j
res[i,j]=1.0
else
res[i,j]=0.0
end
end
end
return res
end
# operator overloading
"""
Base.:*(x::uwreal, y::Matrix{uwreal})
Multiplication operator overloading between uwreal variable and matrix of uwreal
"""
function Base.:*(x::uwreal, y::Matrix{uwreal})
n,m = size(y)
res = Matrix{uwreal}(undef, n,m)
for i in 1:n
for j in 1:m
res[i,j] = x * y[i,j]
if res[i,j] !=0 && res[i,j] !=1
#uwerr(res[i,j])
end
end
end
return res
end
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment