small modification in juobs_linalg and analysis code

gevp tools/tools.jl

@@ -85,9 +85,9 @@ mutable struct  infoMatt
     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))
+        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]
@@ -140,44 +140,43 @@ function comp_mat(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr
     end
     return res
 end
-function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false)
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=3)
     aux_mat = comp_mat(tot11, tot12, tot22)
     y0 = getfield(aux_mat[1], :y0)
     res_en = Array{infoEn}(undef, length(aux_mat))
     if !evec
         for i = 1:length(aux_mat)
-            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
-            evalues = uwgevp_tot(mat_obs, 5)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues = uwgevp_tot(mat_obs, tnot)
             res_en[i] = infoEn(energies(evalues), aux_mat[i])
         end
         return res_en
     else
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
-            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
-            evalues, eigvec = uwgevp_tot(mat_obs, 5, evec=true)
-            println("in comp_energy, i is = ", i)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true, iter=50)
             res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
         return res_en, evec
     end
 end
-function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=5)
+function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.Corr}}, tot13::Array{Array{juobs.Corr}}, tot22::Array{Array{juobs.Corr}}, tot33::Array{Array{juobs.Corr}}, tot23::Array{Array{juobs.Corr}}, evec::Bool=false; tnot::Int64=2)
     aux_mat = comp_mat(tot11, tot12, tot13, tot22, tot23, tot33)
     y0 = getfield(aux_mat[1], :y0)
     res_en = Array{infoEn}(undef, length(aux_mat))
     if !evec
         for i = 1:length(aux_mat)
-            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
-            evalues = uwgevp_tot(mat_obs, tnot)
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
+            evalues = uwgevp_tot(mat_obs, tnot, iter=50)
             res_en[i] = infoEn(energies(evalues), aux_mat[i])
         end
         return res_en
     else
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
-            mat_obs = getfield(aux_mat[i], :mat)[y0:end-1]
+            mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
             evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
             println("in comp_energy, i is = ", i)
             res_en[i] = infoEn(energies(evalues), aux_mat[i])
@@ -197,18 +196,18 @@ mutable struct infoEn
         new(en_val, mu, y0)
     end
 end
-function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, ppcorr::Array{uwreal,1}, mass::uwreal)
-    Rn = [ (uwdot(evec[i][:,1], uwdot(ppcorr[28+i], evec[i][:,1])))^(0.5) * exp(mass*i/2)  for i=1:length(evec)-2]
-    return Rn
+function pseudo_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, mu::Array{Float64})
+    Rn = [exp(mass * (t)/2) * evec[t][1] for t =2:length(evec)]
+    aux = sqrt(2)*sum(mu)  / mass^1.5 
+    return uwdot(aux,Rn)
+end
+function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal)
+    Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
+    return uwdot(Rn , 1/ mass)
 end
 function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
     return sqrt(2)*sum(mu)*mat_elem / mass^1.5
 end
-function vec_mat_elem(evec::Array{Array{T,2} where T,1}, akakcorr::Array{uwreal,1},gevpmat::Array{Array{T,2} where T,1},  mass::uwreal )
-    Rn = [ (uwdot(evec[i][:,2], uwdot(akakcorr[28+i], evec[i][:,2]))^2)^(-0.25) * exp(mass*i/2)  for i=1:length(evec)-2]
-    aux = vec([ uwdot(evec[i][:,2], gevpmat[28+i][:,2]) for i =1:length(evec)-2])
-    return Rn.*aux
-end
 function match_muc(muh, m_lh, m_lh_star, target)
     M = (m_lh .+ 3 .* m_lh_star) ./ 4
 

src/juobs_linalg.jl

@@ -117,40 +117,33 @@ 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)
+    n = length(mat_list)
+    evals = Array{Array{uwreal}}(undef, 0 )
     if !evec
+        println(tnot)
         c_inv = invert(mat_list[tnot])
-        aux_evals = uwgevp.(mat_list, c_tnot_inv=c_inv, iter = iter)
-        evals = get_diag.(aux_evals)
+        for i =1:n
+            aux_evals = uwgevp(mat_list[i], c_inv, iter = iter)
+            push!(evals , get_diag(aux_evals))
+        end
         return evals
     else
-        n = length(mat_list)
-        println(n)
-        evals = Array{Array{uwreal}}(undef, 0 )
         evecs = Array{Matrix}(undef, 0)
-        for i = 6:n-1
-            println("The time t is =", i)
-            #tnot+=1
-            c_inv = invert(mat_list[tnot])
-            println("The time t0 is =", tnot)
-            aux_evals, aux_evec =  uwgevp(mat_list[i], c_tnot_inv=c_inv, iter = iter, evec = evec)
+        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))
-            push!(evecs, aux_evec)
+            ptilda = norm_evec(aux_evec, mat_list[i])
+            m = uwdot(mat_list[i], ptilda)
+            push!(evecs, m)
+            println( "\n t is ", i+3, " tnot is ", i)            
+            #if  2*tnot < i && i<n-5
+            #    tnot+=5
+            #end
         end
-        println("In uwgevp_tot, the time extent is ", length(evals))
-        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_inv::Matrix{uwreal}; iter = 30, evec = false)
 
@@ -160,7 +153,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_inv::Matrix{uwreal}, iter::Int64 = 30, evec::Bool = false)
+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
@@ -191,7 +184,18 @@ function uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
     end
 end
 
-
+function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
+    n=size(evec,1)
+    res_evec=Matrix{uwreal}(undef, n, n)
+    for i =1:n
+        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})
     n = size(a,1)
     res = [a[k, k] for k = 1:n]