improving performance in juobs_linalg by adding specific function for eigenvectors

gevp tools/tools.jl

@@ -138,11 +138,11 @@ function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.C
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
             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])
+            eigvec = uwgevp_tot(mat_obs, delta_t=3)
+            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
-        return res_en, evec
+        return  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=2)
@@ -160,12 +160,12 @@ function comp_energy(tot11::Array{Array{juobs.Corr}}, tot12::Array{Array{juobs.C
         evec = Array{Array{}}(undef, length(aux_mat))
         for i = 1:length(aux_mat)
             mat_obs = getfield(aux_mat[i], :mat)[y0+1:end-1]
-            evalues, eigvec = uwgevp_tot(mat_obs, tnot, evec=true)
+            eigvec = uwgevp_tot(mat_obs, delta_t=3)
             println("in comp_energy, i is = ", i)
-            res_en[i] = infoEn(energies(evalues), aux_mat[i])
+            #res_en[i] = infoEn(energies(evalues), aux_mat[i])
             evec[i] = eigvec
         end
-        return res_en, evec
+        return  evec
     end
 end
 mutable struct infoEn
@@ -180,13 +180,22 @@ mutable struct infoEn
     end
 end
 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)]
+    Rn = [exp(mass * (t)/2) * evec[t][1] for t =1: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)
+function vec_mat_elem(evec::Array{Array{T,2} where T,1}, mass::uwreal, deg::Bool)
+    if !deg
+        Rn = [exp(mass * t/2) * evec[t][1] for t =1:length(evec)]
+        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
+    else
+        Rn = [exp(mass * t/2) * evec[t][4] for t =1:length(evec)]
+        println(za(beta[iens]))
+        println(mass)
+        return uwdot(Rn , za(beta[iens]) * sqrt(2/ mass))
+        #return uwdot(Rn , 1/mass)# * sqrt(2/ mass))
+end
+
 end
 function extract_dec_const(mat_elem::uwreal, mass::uwreal, mu::Array{Float64})
     return sqrt(2)*sum(mu)*mat_elem / mass^1.5

src/juobs_linalg.jl

@@ -120,7 +120,6 @@ function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, eve
     n = length(mat_list)
     evals = Array{Array{uwreal}}(undef, 0 )
     if !evec
-        println(tnot)
         c_inv = invert(mat_list[tnot])
         for i =1:n
             aux_evals = uwgevp(mat_list[i], c_inv, iter = iter)
@@ -136,18 +135,26 @@ function uwgevp_tot(mat_list::Vector{Matrix}, tnot::Int64; iter::Int64 = 30, eve
             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
         return evals, evecs
     end
 end
+
+function uwgevp_tot(mat_list::Vector{Matrix}; iter::Int64=30, delta_t::Int64=3)
+    n = length(mat_list) - delta_t
+    evecs = Array{Matrix}(undef, 0)
+    for i = 2:n
+        c_inv = invert(mat_list[i])
+        usless, temp_evec = uwgevp(mat_list[i+delta_t], c_inv, evec=true, iter=iter)
+        push!(evecs, uwdot(mat_list[i], norm_evec(temp_evec, mat_list[i])))
+    end
+    return evecs
+end
+
 """
 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
+Given two matrices of uwreal C(t) and C(t0)^-1, 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.