implementation of tridiagonal reduction algorithms to improve eigenvalue problem performances

src/juobs.jl

@@ -9,7 +9,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
-export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
+export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction
 export corr_obs, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 

src/juobs_linalg.jl

@@ -126,6 +126,7 @@ function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} }; wp
             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
@@ -173,7 +174,7 @@ corr_upper = [corr_pa[1]]
 matrices = [corr_diag, corr_upper]
 
 ## solve the GEVP
-evals = getall_eigvals(matrices, 5) #where t_0=5
+#evals = getall_eigvals(matrices, 5) #where t_0=5
 
 
 Julia>
@@ -213,23 +214,46 @@ res = uweigvals(a)  ##eigenvalues
 res1 = uweigvals(a,b) ## generalised eigenvalues
 ```
 """
-function uweigvals(a::Matrix{uwreal}; iter = 30)
+function uweigvals(a::Matrix{uwreal}; iter::Int64=5, diag::Bool=true)
     n = size(a,1)
-    for k in 1:iter
-        q_k, r_k = qr(a)
-        a = uwdot(r_k, q_k)
+    if n == 2 
+        return uweigvals_dim_geq3(a, iter=iter, diag=diag) #uweigvals_dim2(a) 
+    else
+        return uweigvals_dim_geq3(a, iter=iter, diag=diag)
     end
-    return get_diag(a)
 end
-function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
+function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 5, diag::Bool=true)
     c = uwdot(invert(b),a)
+    n = size(c,1)
+    if n == 2
+        return uweigvals_dim_geq3(c, iter=iter, diag=diag) #uweigvals_dim2(c) 
+    else
+        return uweigvals_dim_geq3(c, iter=iter, diag=diag)
+    end
+    #for k in 1:iter
+    #    q_k, r_k = qr(c)
+    #    c = uwdot(r_k, q_k)
+    #end
+end
+function uweigvals_dim_geq3(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    n = size(a,1)
+    a = tridiag_reduction(a)
     for k in 1:iter
-        q_k, r_k = qr(c)
-        c = uwdot(r_k, q_k)
+        q_k, r_k = qr(a)
+        a = uwdot(r_k, q_k)
     end
-    return get_diag(c)
+    diag == true ? (return get_diag(a)) : (return a)
+end
+function uweigvals_dim2(a::Matrix{uwreal})
+    res = Vector{uwreal}(undef,2)
+    aux_sum = a[1,1] + a[2,2]
+    delta = ADerrors.sqrt((aux_sum)^2 - 4*(a[1,1]*a[2,2]- a[1,2]*a[2,1]))
+    res[1] = (aux_sum + delta) / 2
+    res[2] = (aux_sum - delta) / 2
+    return res
 end
 
+
 @doc raw"""
     getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30 )
 
@@ -257,7 +281,7 @@ mat_array = get_matrix(diag, upper_diag)
 evecs = getall_eigvecs(mat_array, 5)
 ```
 """
-function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30)
+function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=5)
     n = length(a)-delta_t
     res = Vector{Matrix{uwreal}}(undef, n)
     [res[i] = uweigvecs(a[i+delta_t], a[i]) for i=1:n]
@@ -291,18 +315,22 @@ res = uweigvecs(a)  ##eigenvectors in column
 res1 = uweigvecs(a,b) ## generalised eigenvectors in column 
 ```
 """
-function uweigvecs(a::Matrix{uwreal}; iter = 30)
+function uweigvecs(a::Matrix{uwreal}; iter = 5)
     n = size(a,1)
+    if n > 2
+        a, q = tridiag_reduction(a, q_mat=true)
+    end
     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  u
+    n > 2 ? (return  uwdot(q,u)) : (return u)
 end
-function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5)
     c  = uwdot(invert(b), a)
+    #return uweigvecs(c, iter=iter)
     n = size(c,1)
     u = idty(n)
     for k in 1:iter
@@ -343,11 +371,11 @@ eval, evec = uweigen(a)
 eval1, evec1 = uweigvecs(a,b) 
 ```
 """ 
-function uweigen(a::Matrix{uwreal}; iter = 30)
-    return uweigvals(a, iter=iter), uweigvecs(a, iter=iter)
+function uweigen(a::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    return uweigvals(a, iter=iter,diag=diag), uweigvecs(a, iter=iter)
 end
-function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
-    return uweigvals(a, b, iter=iter), uweigvecs(a, b, iter=iter)
+function uweigen(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 5, diag::Bool=true)
+    return uweigvals(a, b, iter=iter, diag=diag), uweigvecs(a, b, iter=iter)
 end
 
 function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
@@ -379,13 +407,32 @@ Given an input matrix of uwreals, this function returns the inverse of that matr
 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
 
 
+"""
+Performs a Cholesky decomposition of A, which must 
+    be a symmetric and positive definite matrix. The function
+    returns the lower variant triangular matrix, L.
+"""
+function uwcholesky(A::Matrix{uwreal})
+    n = size(A,1)
+    L = fill(uwreal(0), n, n)
+    for i in 1:n 
+        for k in 1:i
+            tmp_sum = sum(L[i,j] * L[k,j] for j in 1:k)
+            if i==k 
+                L[i,k] = sqrt(((A[i,i] - tmp_sum)))
+            else
+                L[i,k] = (1.0 / L[k,k] * (A[i,k]-tmp_sum))
+            end
+        end
+    end
+    return L 
+end
 """
 qr(A::Matrix{uwreal})
 
@@ -441,13 +488,6 @@ function make_householder(a::Vector{uwreal})
     v[1] = uwreal(1.0) 
     H = idty(n)
     H = H -  (2.0 / uwdot(v,v)) * uwdot(reshape(v, :,1), transpose(v)) #this last term is a vector product  
-    for i in 1:n
-        for j in 1:n
-            if H[i,j]!= 0 &&  H[i,j]!=1 
-
-            end
-        end
-    end
     return H
 end
 
@@ -462,28 +502,99 @@ function make_householder(a::Vector{Float64})
     H = H .-  (2.0 / dot(v,v)) * mul!(y, v, transpose(v))
     return H
 end
-
+"""
+This function computes a vector u that generates a Householder reflection H = I - uu* satisfying Ha=nu*e1.
+Accumulating this transformations any matrix can
+be reduced to upper Hessenberg form.
+"""
+function house_gen(a::Vector{uwreal})
+    u = deepcopy(a)
+    nu = uwnorm(a)
+    if nu == 0 
+        u[1] = sqrt(2)
+    end
+    if u[1] != 0 
+        rho = u[1] / sqrt(u[1]*u[1])
+    else
+        rho = 1
+    end
+    u = uwdot(rho/nu, u) 
+    u[1] += 1
+    u =  uwdot(u , 1/sqrt(u[1]))
+    nu = -(rho) * nu
+    return u, nu 
+end
 """
 hess_red(a::Matrix{uwreal})
 
-Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix.
+Given a matrix of uwreal variables, this function returns the Hessenberg form of the matrix by applying householder transformations.
 """
-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]
+function hess_reduce(a::Matrix{uwreal}; q_mat::Bool=false)
+    n = size(a,1)
+    H = deepcopy(a)
+    Q = idty(n)
+    
+    for k in 1:n-2
+        u, H[k+1,k] = house_gen(H[k+1:n,k])
+        Q[k+1:n,k] = u
+
+        v = uwdot(u,H[k+1:n, k+1:n])
+        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
+        H[k+2:n,k] .= 0.0 
+        
+        v = uwdot(H[1:n,k+1:n], u)
+        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
     end
-        res = uwdot(mi_inv, uwdot(res, mi) )
+    
+    if q_mat 
+        Id = idty(n)
+        for k  in  reverse(1:n-2)
+            u = Q[k+1:n, k]
+            v = uwdot(u, Q[k+1:n,k+1:n])
+            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
+            Q[:,k] = Id[:,k]
         end
-    return res
+        return H,Q
+    else
+        return H
+    end
+    
 end
+"""
+
 
+"""
+function tridiag_reduction(a::Matrix{uwreal}; q_mat::Bool=false)
+    n = size(a,1)
+    H = deepcopy(a)
+    Q = idty(n)
+    for k in 1:n-2
+        u, H[k+1,k] = house_gen(H[k+1:n,k])
+        Q[k+1:n,k] = u
+
+        v = uwdot(u,H[k+1:n, k+1:n])
+        H[k+1:n, k+1:n] .= H[k+1:n, k+1:n] .- uwdot(reshape(u,length(u),1),v)
+        H[k+2:n,k] .= 0.0 
+        
+        v = uwdot(H[1:n,k+1:n], u)
+        H[1:n,k+1:n] .= H[1:n,k+1:n] .- uwdot(v,reshape(u,1,length(u)))
+        H[k,k+2:n] .= 0.0 
+    end
+    
+    if q_mat 
+        Id = idty(n)
+        for k  in  reverse(1:n-2)
+            u = Q[k+1:n, k]
+            v = uwdot(u, Q[k+1:n,k+1:n])
+            Q[k+1:n,k+1:n] .= Q[k+1:n,k+1:n] .- uwdot(reshape(u,length(u),1),v)
+            Q[:,k] = Id[:,k]
+        end
+        return H,Q
+    else
+        return H
+    end
+    
+end
 """
 upper_inversion(u::Matrix{uwreal})
 
@@ -556,7 +667,7 @@ function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
     k,p = size(b)
     c = Matrix{uwreal}(undef, n, p)
     if m != k
-        error("Error uwdot")
+        error("Error uwdot: in a*b the size of a is ", size(a), " while the size of b is ", size(b) )
     end
     for i = 1:n
         for j = 1:p
@@ -573,7 +684,7 @@ It returns the norm of a vector of uwreal
 """
 function uwnorm(a::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     norm = sqrt(uwdot(a,a))
-    isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
+    #isnothing(wpm) ? uwerr(norm) : uwerr(norm, wpm)
     return norm 
 end
 
@@ -631,7 +742,6 @@ function idty(n::Int64)
 end
 
 
-# operator overloading
 
 """
 Base.:*(x::uwreal, y::Matrix{uwreal})