introducing new routines for gevp

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
+export get_matrix, uwgevp_tot, energies, uwdot, uweigvals, uweigvecs, uweigen, invert
 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

@@ -11,6 +11,19 @@
 #
 #
 ##################################################################
+using ForwardDiff, LinearAlgebra
+for op in (:eigvals, :eigvecs)
+    @eval function LinearAlgebra.$op(a::Matrix{uwreal})
+
+        function fvec(x::Matrix)
+            return LinearAlgebra.$op(x)
+        end
+        return fvec(value.(a))
+        #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} )
@@ -172,23 +185,51 @@ 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 uwevp(a::Matrix{uwreal}; iter = 30, evec = false)
+function uweigvals(a::Matrix{uwreal}; iter = 30)
     n = size(a,1)
-    if evec
-        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 a, u
-    else
-        for k in 1:iter
-            q_k, r_k = qr(a)
-            a = uwdot(r_k, q_k)
-        end
-        return a
+    for k in 1:iter
+        q_k, r_k = qr(a)
+        a = uwdot(r_k, q_k)
+    end
+    return a#get_diag(a)
+end
+function uweigvals(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+    c  = uwdot(invert(b), a)
+    n = size(c,1)
+    for k in 1:iter
+        q_k, r_k = qr(c)
+        c = uwdot(r_k, q_k)
     end
+    return c#get_diag(c)
+end
+
+
+function uweigvecs(a::Matrix{uwreal}; iter = 30)
+    n = size(a,1)
+    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
+end
+function uweigvecs(a::Matrix{uwreal}, b::Matrix{uwreal}; iter = 30)
+    c  = uwdot(invert(b), a)
+    n = size(c,1)
+    u = idty(n)
+    for k in 1:iter
+        q_k, r_k = qr(c)
+        c = uwdot(r_k, q_k)
+        u = uwdot(u, q_k)
+    end
+    return  u
+end
+function uweigen(a::Matrix{uwreal}; iter = 30)
+    return uweigvals(a, iter=iter), 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)
 end
 
 function norm_evec(evec::Matrix{uwreal}, ctnot::Matrix{uwreal})
@@ -414,7 +455,7 @@ It returns the norm of a vector of uwreal
 function uwnorm(a::Vector{uwreal})
     norm = sqrt(uwdot(a,a))
     uwerr(norm)
-    return norm
+    return norm 
 end
 
 """