First working case

Author: albertomercurio

src/qobj/eigsolve.jl

@@ -144,11 +144,44 @@ function _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, k, m, β, sorted_vals, sort
     ordschur!(F, select)
 
     copyto!(Hₘ, F.T)
-    Tₘ = Hₘ
     copyto!(Uₘ, F.Z)
     mul!(f, Uₘᵥ, β)
 
-    return Tₘ, Uₘ
+    return nothing
+end
+
+# Pure Julia implementation of computing right eigenvectors from Schur form
+# Instead of using LAPACK.trevc!('R', 'A', select, Tₘ)
+function _schur_right_eigenvectors(Tₘ, k)
+    n = size(Tₘ, 1)
+    vecs = zeros(eltype(Tₘ), n, k)
+    k == 0 && return vecs
+
+    value_tol(λ) = eps(typeof(abs(λ))) * max(one(typeof(abs(λ))), abs(λ))
+
+    @inbounds for col in 1:k
+        vec = view(vecs, :, col)
+        fill!(vec, zero(eltype(Tₘ)))
+        vec[col] = one(eltype(Tₘ))
+        λ = Tₘ[col, col]
+
+        for row in (col-1):-1:1
+            acc = zero(eltype(Tₘ))
+            for inner in (row + 1):col
+                acc += Tₘ[row, inner] * vec[inner]
+            end
+            denom = Tₘ[row, row] - λ
+            if abs(denom) <= value_tol(λ)
+                vec[row] = zero(eltype(Tₘ))
+            else
+                vec[row] = -acc / denom
+            end
+        end
+
+        LinearAlgebra.normalize!(vec)
+    end
+
+    return vecs
 end
 
 function _eigsolve(
@@ -196,12 +229,13 @@ function _eigsolve(
     V₁ₖ = view(V, :, 1:k)
     Vₖ₊₁ = view(V, :, k + 1)
     Hₖ₊₁₁ₖ = view(H, k + 1, 1:k)
+    cache0₁ₖ = view(cache0, :, 1:k)
     cache1₁ₖ = view(cache1, :, 1:k)
     cache2₁ₖ = view(cache2, 1:k)
 
     M = typeof(cache0)
 
-    Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, k, m, β, sorted_vals, sortby, rev)
+    _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, k, m, β, sorted_vals, sortby, rev)
 
     numops = m
     iter = 0
@@ -227,20 +261,18 @@ function _eigsolve(
 
         # println( A * Vₘ ≈ Vₘ * M(Hₘ) + qₘ * M(transpose(βeₘ)) )     # SHOULD BE TRUE
 
-        Tₘ, Uₘ = _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, k, m, β, sorted_vals, sortby, rev)
+        _update_schur_eigs!(Hₘ, Uₘ, Uₘᵥ, f, k, m, β, sorted_vals, sortby, rev)
 
         numops += m - k - 1
         iter += 1
     end
 
+    Tₘ = Hₘ
     vals = diag(view(Tₘ, 1:k, 1:k))
-    select = Vector{Int}(undef, 0)
-    VR = LAPACK.trevc!('R', 'A', select, Tₘ)
-    @inbounds for i in 1:size(VR, 2)
-        normalize!(view(VR, :, i))
-    end
-    mul!(cache1, Vₘ, M(Uₘ * VR))
-    vecs = cache1[:, 1:k]
+    VR = _schur_right_eigenvectors(Tₘ, k)
+    mul!(cache0₁ₖ, Uₘ, VR)
+    mul!(cache1₁ₖ, Vₘ, cache0₁ₖ)
+    vecs = copy(cache1₁ₖ)
     settings.auto_tidyup && tidyup!(vecs)
 
     return EigsolveResult(vals, vecs, type, dimensions, iter, numops, (iter < maxiter))

test_env/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"

test_env/my_test.jl

@@ -0,0 +1,55 @@
+using Revise
+using LinearAlgebra
+using SparseArrays
+using QuantumToolbox
+# using GenericLinearAlgebra
+using GenericSchur
+using BenchmarkTools
+
+# include("arnoldi_new.jl")
+
+# %%
+
+# T = ComplexF64
+T = Complex{BigFloat}
+
+# A = rand(T, 100, 100)
+A = sprand(T, 100, 100, 0.01) + Diagonal(rand(T, 100))
+x = rand(T, 100)
+
+nev = 5
+krylov_dim = 15
+tol = 1e-8
+
+# values0 = eigvals(Array(A); sortby=abs2)[end-nev+1:end] |> reverse
+res0 = eigen(Array{ComplexF64}(A); sortby=abs2)
+values0 = res0.values[end-nev+1:end]
+vectors0 = res0.vectors[:, end-nev+1:end]
+
+res = eigsolve(A; v0=x, eigvals=nev, krylovdim=krylov_dim, eigstol=tol);
+res.iter
+
+values0
+res.values
+
+idxs1 = sortperm(res.values, by=abs2)
+idxs2 = sortperm(values0, by=abs2)
+
+res.values[idxs1] - values0[idxs2]
+
+# res.vectors[:, idxs1] - vectors0[:, idxs2]
+
+values_2 = [v' * A * v for v in eachcol(res.vectors[:, idxs1])]
+
+values_2 - values0[idxs2]
+
+map(1:nev) do i
+    v0 = vectors0[:, idxs2[i]] |> copy
+    v = res.vectors[:, idxs1[i]] |> copy
+
+    abs(dot(v0, v))
+end
+
+# %%
+
+@benchmark eigsolve($A; v0=$x, eigvals=$nev, krylovdim=$krylov_dim, eigstol=$tol)