From 95ccd1105771953cdbd31b447cf393fdf1b3b5d1 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Fri, 27 Dec 2024 18:17:41 -0500
Subject: [PATCH] Gram-Schmidt implemented

---
 helper.jl    | 25 ++++++++++++++++++++++++-
 test/misc.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 68 insertions(+), 1 deletion(-)
 create mode 100644 test/misc.jl

diff --git a/helper.jl b/helper.jl
index e48b171..6f52421 100644
--- a/helper.jl
+++ b/helper.jl
@@ -47,4 +47,27 @@ function exportCSV(file::String, data, labels=nothing)
             writedlm(f, hcat(reim(d)..., l))
         end
     end
-end
\ No newline at end of file
+end
+
+"c-orthonormalization via Gram-Schmidt"
+function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
+    function proj(ui, vi)
+        num = sum(vecs[ui] .* ws .* vecs[vi])
+        den = sum(vecs[ui] .* ws .* vecs[ui])
+        return (num / den) .* vecs[ui]
+    end
+    mask = ones(Bool, size(vecs)) # vectors to retain
+    for i in eachindex(vecs)
+        for j in findall(mask[1 : (i - 1)])
+            vecs[i] -= proj(j, i)
+        end
+        norm = sum(vecs[i] .* ws .* vecs[i])
+        verbose && println("Mag. of basis vector $i = $(abs(norm))")
+        if abs(norm) > precision
+            vecs[i] = vecs[i] ./ sqrt(norm)
+        else
+            mask[i] = false
+        end
+    end
+    return vecs[mask]
+end
diff --git a/test/misc.jl b/test/misc.jl
new file mode 100644
index 0000000..c7126d8
--- /dev/null
+++ b/test/misc.jl
@@ -0,0 +1,44 @@
+using SparseArrays, LinearAlgebra
+include("../helper.jl")
+
+#### gram_schmidt ####
+ 
+n = 64
+N = 8
+
+vecs = [rand(n) + 1im .* rand(n) for _ in 1:N]
+
+ws = rand(n)
+ws_mat = spdiagm(ws)
+
+basis = hcat(vecs...)
+
+H = rand(n, n)
+H += transpose(H) # complex symmetric
+
+H_EC = transpose(basis) * ws_mat * H * ws_mat * basis
+N_EC = transpose(basis) * ws_mat * basis
+
+evals = eigvals(H_EC, N_EC)
+println("Eigenvalues with GEVP:")
+display(evals)
+
+push!(vecs, vecs[end]) # duplicate last basis vector to test orthogonalization
+ortho_vecs = gram_schmidt(vecs, ws)
+@assert length(ortho_vecs) == N "Basis dimensionality mismatch after Gram-Schmidt"
+ortho_basis = hcat(ortho_vecs...)
+
+N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis
+println("Norm matrix after Gram-Schmidt:")
+display(round.(N_ortho; digits=2)) # should be ≈I
+@assert N_ortho ≈ I(N) "Gram-Schmidt did not yield an orthogonal basis"
+
+H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
+
+evals_ortho = eigvals(H_EC_ortho)
+println("Eigenvalues after Gram-Schmidt:")
+display(evals_ortho)
+
+@assert evals ≈ evals_ortho "Gram-Schmidt did not preserve the eigenvalues"
+
+######################
\ No newline at end of file