Repository navigation
Conversation
GMRES-DR was added to matlab/src/solvers but never listed alongside the other three solvers. Note it as MATLAB-only until a Julia port exists.
Port gmres_dr and its qrupdate_gs incremental-QR utility to Julia, following the existing gmres_e.jl/harmonic_ritz_vectors.jl idioms. Unlike the MATLAB implementation's hybrid eigs/Schur harmonic Ritz strategy (needed to work around ARPACK's eigs being unreliable on small/ill-conditioned matrices), the Julia port uses the dense generalized eigenproblem eigen(F, G) uniformly -- consistent with harmonic_ritz_vectors.jl and migration_plan.md's existing rationale for avoiding Arpack.jl. Complex-conjugate harmonic Ritz pairs are handled inline (both real and imaginary parts recycled) rather than via a separate Schur fallback path. Caught and fixed one bug during porting: F*y = lambda*G*y is a reciprocal formulation (Morgan 2000, eq. 2.10), so the harmonic Ritz values to deflate correspond to the LARGEST-magnitude eigenvalues of this pencil, not the smallest. Verified against MATLAB: sherman1 (m=27, k=3) converges in 72 cycles in Julia vs 94 in MATLAB to the same 1e-12 tolerance; embree3 (m=2, k=1) reproduces MATLAB's documented stall exactly (relres 0.4629, same x). Also add test_qrupdate_gs.jl and test_gmres_dr.jl (mirroring the MATLAB suites), register gmres_dr in KrySBAS.jl, and document it in the Julia docs and top-level README. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CI's julia_style check fails on the whole julia/ directory if any single file is non-compliant. benchmark.jl, diagnose_lgmres.jl, profile_solvers.jl, lgmres.jl, and test_poisson.jl had stray trailing semicolons predating this branch's GMRES-DR work; strip them so the format check passes again. No semantic change (semicolons are inert line punctuation in Julia scripts, not REPL output suppression).
jccf19
left a comment
jccf19
left a comment
There was a problem hiding this comment.
I have added three questions regarding the numerical consistency and robustness of the GMRES-DR implementation.
| # Dispatch: m == n → full unrestarted GMRES | ||
| if m == n | ||
| t0 = time() | ||
| x, stats = Krylov.gmres(A, b; memory = n) |
There was a problem hiding this comment.
Question — Initial guess consistency
It seems that x_initial is not propagated to the fallback Krylov.gmres call. However, res0 below is computed using x_initial. Could you please verify whether x_initial is being properly handled in this case?
| kdvec = Int[] | ||
|
|
||
| V = zeros(T, n, m + 1) | ||
| V[:, 1] = r0 / beta |
There was a problem hiding this comment.
Question — Zero initial residual
If x_initial already satisfies the linear system, then r0 = b - A*x_initial and consequently beta = norm(r0) is zero. In that case, this normalization would involve a division by zero. Could you please verify whether the case beta == 0 is handled before starting the Arnoldi process?
| w .-= H[i, j] .* view(V, :, i) | ||
| end | ||
| H[j+1, j] = norm(w) | ||
| V[:, j+1] = w / H[j+1, j] |
There was a problem hiding this comment.
Question — Arnoldi breakdown handling
Since H[j+1, j] = norm(w), this value may become zero (or numerically very small) in the case of an Arnoldi happy breakdown. Then the normalization w / H[j+1, j] would involve a division by zero. Could you please verify whether this situation is explicitly handled before normalizing the next Arnoldi vector?
This PR introduces the Julia version of GMRES-DR with appropriate unit and functional tests.