Add retrying normalization if necessary.

See JuliaLinearAlgebra#25.

Project.toml

@@ -4,7 +4,9 @@ version = "0.4.3"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 Adapt = "2.1, 3.0"

src/svd.jl

@@ -1,3 +1,5 @@
+using KrylovKit
+using Random: rand!
 function _biLanczosIterations!(A, stepsize, αs, βs, U, V, μs, νs, maxνs, maxμs, τ, reorth_in, tolreorth, debug)
 
     m, n = size(A)
@@ -50,6 +52,20 @@ function _biLanczosIterations!(A, stepsize, αs, βs, U, V, μs, νs, maxνs, ma
                 nReorthVecs += 1
             end
             α = norm(v)
+
+            if α < τ * maximum(size(A)) * eps()  # orthogonalization failed, see https://github.com/poulson/PROPACK/blob/2465f89d5b1fba56de71c3e69e27d017c3dc2295/double/dlanbpro.F#L384
+                @warn "Restart orthogonalization"
+                b = V[1:(j-1)]
+                B = KrylovKit.OrthonormalBasis(b ./ norm.(b))
+                for _ in 1:3
+                    v = rand!(v) * 2 .- 1
+                    KrylovKit.orthonormalize!(v, B, KrylovKit.ModifiedGramSchmidt())
+                    α = norm(v)
+                    if !(α < τ * maximum(size(A)) * eps())
+                        break
+                    end
+                end
+            end
         end
 
         ## update the result vectors