diff --git a/calculations/B2R_Berggren_poles.jl b/calculations/B2R_Berggren_poles.jl
index 3ab147c..76c1a43 100644
--- a/calculations/B2R_Berggren_poles.jl
+++ b/calculations/B2R_Berggren_poles.jl
@@ -33,7 +33,7 @@ for (j, c) in enumerate(training_points)
 end
 
 EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
-N_EC = transpose(EC_basis) * EC_basis
+EC_basis = gram_schmidt!(EC_basis)
 
 extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
 
@@ -47,7 +47,7 @@ for (j, c) in enumerate(extrapolate_points)
     H = get_H_matrix(V_system(c), p, w)
     H_berg = transpose(berg_basis_w) * H * berg_basis
     H_EC = transpose(EC_basis) * H_berg * EC_basis
-    evals = eigvals(H_EC, N_EC)
+    evals = eigvals(H_EC)
     i = argmin(abs.(evals .- exact_E[j]))
     extrapolate_E[j] = evals[i]
 end
diff --git a/calculations/B2R_comparison.jl b/calculations/B2R_comparison.jl
index 7b5d3c6..2459e07 100644
--- a/calculations/B2R_comparison.jl
+++ b/calculations/B2R_comparison.jl
@@ -1,4 +1,5 @@
 using Plots
+include("../helper.jl")
 include("../p_space.jl")
 
 berggren_mesh = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
@@ -23,8 +24,8 @@ for mesh in (berggren_mesh, csm_mesh)
     end
 
     EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
+    EC_basis = gram_schmidt!(EC_basis, w)
     EC_basis_w = EC_basis .* w
-    N_EC = transpose(EC_basis_w) * EC_basis
 
     extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
 
@@ -36,7 +37,7 @@ for mesh in (berggren_mesh, csm_mesh)
 
         H = get_H_matrix(V_system(c), p, w)
         H_EC = transpose(EC_basis_w) * H * EC_basis
-        evals = eigvals(H_EC, N_EC)
+        evals = eigvals(H_EC)
         i = argmin(abs.(evals .- exact_E[j]))
         extrapolate_E[j] = evals[i]
     end
diff --git a/calculations/R2R_Berggren.jl b/calculations/R2R_Berggren.jl
index 8bc17b2..f5656f4 100644
--- a/calculations/R2R_Berggren.jl
+++ b/calculations/R2R_Berggren.jl
@@ -1,4 +1,5 @@
 using Plots
+include("../helper.jl")
 include("../p_space.jl")
 
 # contour
@@ -27,15 +28,15 @@ ref_E = 0.2 - 0.1im
 exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
 extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
 
+EC_basis = gram_schmidt!(EC_basis, w)
 EC_basis_w = EC_basis .* w
-N_EC = transpose(EC_basis_w) * EC_basis
 
 for (j, c) in enumerate(extrapolate_points)
     exact_E[j] = quick_pole_E(V_system(c))
 
     H = get_H_matrix(V_system(c), p, w)
     H_EC = transpose(EC_basis_w) * H * EC_basis
-    evals = eigvals(H_EC, N_EC)
+    evals = eigvals(H_EC)
     i = argmin(abs.(evals .- ref_E))
     global ref_E = evals[i]
     extrapolate_E[j] = evals[i]
diff --git a/calculations/R2R_Berggren_poles.jl b/calculations/R2R_Berggren_poles.jl
index 3586afe..1d68787 100644
--- a/calculations/R2R_Berggren_poles.jl
+++ b/calculations/R2R_Berggren_poles.jl
@@ -33,7 +33,7 @@ for (j, c) in enumerate(training_points)
     EC_basis[:, j] = evecs[:, i]
 end
 
-N_EC = transpose(EC_basis) * EC_basis
+EC_basis = gram_schmidt!(EC_basis)
 
 extrapolate_points = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8)
 
@@ -47,7 +47,7 @@ for (j, c) in enumerate(extrapolate_points)
     H = get_H_matrix(V_system(c), p, w)
     H_berg = transpose(berg_basis_w) * H * berg_basis
     H_EC = transpose(EC_basis) * H_berg * EC_basis
-    evals = eigvals(H_EC, N_EC)
+    evals = eigvals(H_EC)
     i = argmin(abs.(evals .- exact_E[j]))
     extrapolate_E[j] = evals[i]
 end
diff --git a/calculations/XZ_technique.jl b/calculations/XZ_technique.jl
index 7ae8efe..98389e3 100644
--- a/calculations/XZ_technique.jl
+++ b/calculations/XZ_technique.jl
@@ -1,4 +1,5 @@
 using LinearAlgebra, Plots
+include("../helper.jl")
 include("../ho_basis.jl")
 include("../p_space.jl")
 
@@ -36,14 +37,14 @@ for (j, c) in enumerate(train_cs)
     EC_basis[:, j] = evecs[:, i]
 end
 
-N_EC = transpose(EC_basis) * EC_basis
+EC_basis = gram_schmidt!(EC_basis)
 
 for (j, c) in enumerate(target_cs)
     exact_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.5)
 
     H = T + c .* V
     H_EC = transpose(EC_basis) * H * EC_basis
-    evals = eigvals(H_EC, N_EC)
+    evals = eigvals(H_EC)
     i = argmin(abs.(evals .- exact_E[j]))
     extrapolate_E[j] = evals[i]
 end
diff --git a/helper.jl b/helper.jl
index 359948d..69973c6 100644
--- a/helper.jl
+++ b/helper.jl
@@ -51,7 +51,7 @@ end
 
 "In-place c-orthogonalization via (modified) Gram-Schmidt. Only significant vectors are returned (c-normalized).
 The number of significant vectors to return are determined by the original singular values (compared to the threshold)."
-function gram_schmidt!(vecs, ws=ones(length(vecs[1])), threshold=1e-5; verbose=false)
+function gram_schmidt!(vecs::Vector{Vector{T}}, ws=ones(length(vecs[1])), threshold=1e-5; verbose=false) where T<: Number
     c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
     norm(i) = c_product(i, i)
     proj(i, j) = (c_product(i, j) / norm(i)) .* vecs[i] # component of vec[i] in vec[j]
@@ -86,6 +86,9 @@ function gram_schmidt!(vecs, ws=ones(length(vecs[1])), threshold=1e-5; verbose=f
     return vecs[selected_vecs_i]
 end
 
+"Same as above but the basis is provided as a matrix instead of an array of vectors"
+gram_schmidt!(vecs::Matrix{T}, ws=ones(size(vecs, 1)), threshold=1e-5; verbose=false) where T<: Number = hcat(gram_schmidt!([vecs[:, i] for i in axes(vecs, 2)], ws, threshold; verbose=verbose)...)
+
 "Rank of a basis set (pre-normalized) under c-product"
 c_rank(vecs, ws, threshold=1e-8) = count(c_singular_values(vecs, ws) .> threshold)