improving juobs_linalg with new functions for computing eigenvalues and eigenvectors

src/juobs.jl

@@ -8,7 +8,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md
-export get_matrix, uwgevp_tot, energies, uwdot, uweigvals, uweigvecs, uweigen, invert
+export get_matrix, uwgevp_tot, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs
 export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc, comp_t0
 

src/juobs_linalg.jl

-##################################################################
-# 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
-#
-#
-##################################################################
-using ForwardDiff, LinearAlgebra
+#=using ForwardDiff, LinearAlgebra
 for op in (:eigvals, :eigvecs)
     @eval function LinearAlgebra.$op(a::Matrix{uwreal})
 
@@ -22,7 +9,7 @@ for op in (:eigvals, :eigvecs)
         #return uwreal(LinearAlgebra.$op(a.mean), a.prop, ForwardDiff.derivative($op, a.mean)*a.der)
     end
 end
-
+=#
 
 @doc raw"""
 get_matrix(corr_diag::Vector{Array}, corr_upper::Vector{Array} )
@@ -64,6 +51,7 @@ 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
@@ -184,6 +172,12 @@ 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 )
+    n = length(a)
+    res = Vector{Vector{uwreal}}(undef, n)
+    [res[i] = uweigvals(a[i], a[t0]) for i=1:n]
+    return res
+end
 function uweigvals(a::Matrix{uwreal}; iter = 30)
     n = size(a,1)
     for k in 1:iter
@@ -193,8 +187,7 @@ function uweigvals(a::Matrix{uwreal}; iter = 30)
     return get_diag(a)
 end
 function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal};  iter = 30)
-    c  = uwdot(invert(b), a)
-    n = size(c,1)
+    c = uwdot(invert(b),a)
     for k in 1:iter
         q_k, r_k = qr(c)
         c = uwdot(r_k, q_k)
@@ -202,6 +195,13 @@ function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
     return get_diag(c)
 end
 
+function getall_eigvecs(a::Vector{Matrix}, delta_t; iter=30)
+    n = length(a)-delta_t
+    res = Vector{Matrix{uwreal}}(undef, n)
+    [res[i] = uweigvecs(a[i+delta_t], a[i]) for i=1:n]
+    return res 
+
+end
 function uweigvecs(a::Matrix{uwreal}; iter = 30)
     n = size(a,1)
     u = idty(n)