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"

######################