From 56eb7ea2a639c9552215a71aa31c3d9f92568358 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 30 Sep 2026 15:50:28 +0200 Subject: [PATCH 1/3] Size the `BlockLanczos` thick restart by R_size instead of the initial block size --- src/eigsolve/blocklanczos.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/eigsolve/blocklanczos.jl b/src/eigsolve/blocklanczos.jl index 769b2313..746134e1 100644 --- a/src/eigsolve/blocklanczos.jl +++ b/src/eigsolve/blocklanczos.jl @@ -67,21 +67,23 @@ function eigsolve( numops += fact.R_size else # Shrink and restart following the shrinking method of `Lanczos`. numiter >= maxiter && break - keep = max(div(3 * krylovdim + 2 * converged, 5 * bs), 1) * bs - H = zeros(eltype(fact.H), keep + bs, keep) - # The last bs rows of U contribute to calculate errors of Ritz values. + # Use the current residual block size `bs_R`, which can be smaller than `bs`: + # only the last `bs_R` rows of U couple to the residual. + keep = max(div(3 * krylovdim + 2 * converged, 5 * bs_R), 1) * bs_R + H = zeros(eltype(fact.H), keep + bs_R, keep) + # The last bs_R rows of U contribute to calculate errors of Ritz values. @inbounds for j in 1:keep H[j, j] = D[j] - H[(keep + 1):end, j] = view(U, (K - bs + 1):K, j) + H[(keep + 1):end, j] = view(U, (K - bs_R + 1):K, j) end # Turn diagonal matrix D into a block tridiagonal matrix, and make sure # The residual of krylov subspace keeps the form of [0,..,0,R] @inbounds for j in keep:-1:1 - h, ν = householder(H, j + bs, 1:j, j) - H[j + bs, j] = ν - H[j + bs, 1:(j - 1)] .= zero(eltype(H)) + h, ν = householder(H, j + bs_R, 1:j, j) + H[j + bs_R, j] = ν + H[j + bs_R, 1:(j - 1)] .= zero(eltype(H)) lmul!(h, H) - rmul!(view(H, 1:(j + bs - 1), :), h') + rmul!(view(H, 1:(j + bs_R - 1), :), h') rmul!(U, h') end # Transform the basis and update the residual and update the BTD. @@ -92,7 +94,7 @@ function eigsolve( basistransform!(B, view(U, :, 1:keep)) R_new = OrthonormalBasis(fact.R.vec[1:bs_R]) - view_H = view(H, (keep + bs - bs_R + 1):(keep + bs), (keep - bs_R + 1):keep) + view_H = view(H, (keep + 1):(keep + bs_R), (keep - bs_R + 1):keep) basistransform!(R_new, view_H) fact.R.vec[1:bs_R] = R_new[1:bs_R] From 17acdcb1581ab607eab9603d38036515aef19611 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 30 Sep 2026 20:30:00 +0200 Subject: [PATCH 2/3] Don't let `BlockLanczos` restarts keep fewer vectors than have converged --- src/eigsolve/blocklanczos.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/eigsolve/blocklanczos.jl b/src/eigsolve/blocklanczos.jl index 746134e1..f8c9f124 100644 --- a/src/eigsolve/blocklanczos.jl +++ b/src/eigsolve/blocklanczos.jl @@ -70,6 +70,9 @@ function eigsolve( # Use the current residual block size `bs_R`, which can be smaller than `bs`: # only the last `bs_R` rows of U couple to the residual. keep = max(div(3 * krylovdim + 2 * converged, 5 * bs_R), 1) * bs_R + # Rounding down to a multiple of the block size can drop converged Ritz vectors; + # then use the unrounded `Lanczos` value, which satisfies `converged <= keep < krylovdim`. + keep < converged && (keep = div(3 * krylovdim + 2 * converged, 5)) H = zeros(eltype(fact.H), keep + bs_R, keep) # The last bs_R rows of U contribute to calculate errors of Ritz values. @inbounds for j in 1:keep From 7db5f8a8c9eb4aa1e4a979b3d32cdb6d0f533632 Mon Sep 17 00:00:00 2001 From: Jutho Date: Tue, 6 Oct 2026 00:10:52 +0200 Subject: [PATCH 3/3] Remove comment --- src/eigsolve/blocklanczos.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/eigsolve/blocklanczos.jl b/src/eigsolve/blocklanczos.jl index f8c9f124..88742e8e 100644 --- a/src/eigsolve/blocklanczos.jl +++ b/src/eigsolve/blocklanczos.jl @@ -67,8 +67,6 @@ function eigsolve( numops += fact.R_size else # Shrink and restart following the shrinking method of `Lanczos`. numiter >= maxiter && break - # Use the current residual block size `bs_R`, which can be smaller than `bs`: - # only the last `bs_R` rows of U couple to the residual. keep = max(div(3 * krylovdim + 2 * converged, 5 * bs_R), 1) * bs_R # Rounding down to a multiple of the block size can drop converged Ritz vectors; # then use the unrounded `Lanczos` value, which satisfies `converged <= keep < krylovdim`.