bug fixing in eigenvector computation (juobs_linalg), exporting uwdot

gevp tools/tools.jl

@@ -82,24 +82,34 @@ mutable struct  infoMatt
     mat::Vector{Matrix}
     mu::Vector{Float64}
     y0::Int64
-    function infoMatt(a11::juobs.Corr, a12::Array{juobs.Corr}, a22::Array{juobs.Corr})
-        mu = getfield(a11, :mu)
-        y0 = Int64(getfield(a11, :y0))
-        diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
-        subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
+    function infoMatt(a11::Array{juobs.Corr}, a12::Array{juobs.Corr}, a22::Array{juobs.Corr})
+        mu = getfield(a11[1], :mu)
+        y0 = Int64(getfield(a11[1], :y0))
+        elem11 = sum(a11[i].obs for i =1:length(a11))
+        elem12 = sum(a12[i].obs for i =1:length(a12))
+        elem22 = sum(a22[i].obs for i =1:length(a22))
+        diag = [elem11, elem22]
+        subdiag =[ elem12 ]
+        #diag = [a11.obs, a22[1].obs.+a22[2].obs.+a22[3].obs]
+        #subdiag = [a12[1].obs.+a12[2].obs.+a12[3].obs]
         mat = get_matrix(diag, subdiag)
         return new(mat, mu, y0)
     end
 end
-function comp_mat(tot11::Array{juobs.Corr}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
-    mu = getfield.(tot11,:mu)
+function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
+    mu = getfield.(tot11[1],:mu)
     res = Array{infoMatt}(undef, length(mu))
     for i = 1:length(mu)
-        res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
+        println("tot11 size is ", size(tot11,1), " the type is ", typeof(tot11))
+        a11 = [ tot11[j][i] for j = 1:size(tot11,1) ]
+        a12 = [ tot12[j][i] for j = 1:size(tot12,1) ]
+        a22 = [ tot22[j][i] for j = 1:size(tot22,1) ]
+        res[i] = infoMatt(a11, a12, a22)
+        #res[i]=  infoMatt(tot11[i], [tot12[1][i],tot12[2][i],tot12[3][i]], [tot22[1][i],tot22[2][i],tot22[3][i]])
     end
     return res
 end
-function comp_energy(tot11::Array{juobs.Corr}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}})
     aux_mat = comp_mat(tot11, tot12, tot22)
     y0 = getfield(aux_mat[1], :y0)
     res_en = Array{infoEn}(undef, length(aux_mat))

src/juobs.jl

@@ -8,7 +8,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1
-export get_matrix, uwgevp_tot, energies
+export get_matrix, uwgevp_tot, energies, uwdot
 export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc 
 

src/juobs_linalg.jl

@@ -85,7 +85,7 @@ Given a vector where each entry evals[t] is a uwreal array of eigenvalues, this
 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{Vector{uwreal}})
+function energies(evals::Union{Vector{Vector{uwreal}},Array{Array{uwreal}} })
     time = length(evals) 
     n = length(evals[1])
     eff_en = Array{Array{uwreal}}(undef, n)
@@ -123,16 +123,30 @@ function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, eve
         evals = get_diag.(aux_evals)
         return evals
     else
-        aux_evals, aux_evec =  uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter, evec= true)
+        n = length(mat_list)
+        evals = Array{Array{uwreal}}(undef, n )
+        evecs = Array{Matrix}(undef, n)
+        for i = 1:n
+            aux_evals, aux_evec =  uwgevp(mat_list[i], c_tnot_inv=c_inv, iter = iter, evec = evec)
+            evals[i] =  get_diag(aux_evals)
+            evecs[i] = aux_evec
+        end
+        return evals, evecs
+        #=
+        (aux_evals, aux_evec) =  uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter, evec= true)
+        println(typeof(aux_evals))
+        println(typeof(aux_evec))
+
         evals = get_diag.(aux_evals)
         evecs = aux_evec
         return evals, evecs
+        =#
     end
 end
 
 
 """
-uwgevp(c_t:: Matrix{uwreal}, c_tnot::Matrix{uwreal}; iter = 30, evec = false)
+uwgevp(c_t:: Matrix{uwreal}, c_tnot_inv::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.
@@ -354,6 +368,8 @@ of uwreal data type.
 
 """
 uwdot(a::Vector{uwreal}, b::Vector{uwreal}) = sum(a .* b)
+uwdot(a::Vector{uwreal}, b::uwreal) = [a[i] * b for i=1:length(a)]
+uwdot(a::uwreal, b::Vector{uwreal}) = uwdot(b,a) 
 
 function uwdot(a::Matrix{uwreal}, b::Matrix{uwreal})
     n,m = size(a)