Merge branch 'JuliaSmoothOptimizers:main' into Preconditioning

.buildkite/pipeline.yml

@@ -98,3 +98,17 @@ steps:
         Pkg.instantiate()
         include("test/cpu/component_arrays.jl")'
     timeout_in_minutes: 30
+
+  - label: "CPUs -- LinearOperators.jl"
+    plugins:
+      - JuliaCI/julia#v1:
+          version: "1.10"
+    agents:
+      queue: "juliaecosystem"
+    command: |
+      julia --color=yes --project -e '
+        using Pkg
+        Pkg.add("LinearOperators")
+        Pkg.instantiate()
+        include("test/cpu/linear_operators.jl")'
+    timeout_in_minutes: 30

docs/src/preconditioners.md

@@ -30,12 +30,13 @@ Krylov.jl supports both approaches thanks to the argument `ldiv` of the Krylov s
 ## How to use preconditioners in Krylov.jl?
 
 !!! info
-    - A preconditioner only need support the operation `mul!(y, P⁻¹, x)` when `ldiv=false` or `ldiv!(y, P, x)` when `ldiv=true` to be used in Krylov.jl.
+    - A preconditioner only needs to support the operation `mul!(y, P⁻¹, x)` when `ldiv=false` or `ldiv!(y, P, x)` when `ldiv=true` to be used in Krylov.jl.
+    - Additional support for `adjoint` with preconditioners is required in the methods [`BILQ`](@ref bilq) and [`QMR`](@ref qmr).
     - The default value of a preconditioner in Krylov.jl is the identity operator `I`.
 
 ### Square non-Hermitian linear systems
 
-Methods concerned: [`CGS`](@ref cgs), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`BLOCK-GMRES`](@ref block_gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom).
+Methods concerned: [`CGS`](@ref cgs), [`BILQ`](@ref bilq), [`QMR`](@ref qmr), [`BiCGSTAB`](@ref bicgstab), [`DQGMRES`](@ref dqgmres), [`GMRES`](@ref gmres), [`BLOCK-GMRES`](@ref block_gmres), [`FGMRES`](@ref fgmres), [`DIOM`](@ref diom) and [`FOM`](@ref fom).
 
 A Krylov method dedicated to non-Hermitian linear systems allows the three variants of preconditioning.
 

src/bilq.jl

@@ -15,8 +15,9 @@ export bilq, bilq!
 """
     (x, stats) = bilq(A, b::AbstractVector{FC};
                       c::AbstractVector{FC}=b, transfer_to_bicg::Bool=true,
-                      atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
-                      timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
+                      M=I, N=I, ldiv::Bool=false, atol::T=√eps(T),
+                      rtol::T=√eps(T), itmax::Int=0, timemax::Float64=Inf,
+                      verbose::Int=0, history::Bool=false,
                       callback=solver->false, iostream::IO=kstdout)
 
 `T` is an `AbstractFloat` such as `Float32`, `Float64` or `BigFloat`.
@@ -30,6 +31,7 @@ Solve the square linear system Ax = b of size n using BiLQ.
 BiLQ is based on the Lanczos biorthogonalization process and requires two initial vectors `b` and `c`.
 The relation `bᴴc ≠ 0` must be satisfied and by default `c = b`.
 When `A` is Hermitian and `b = c`, BiLQ is equivalent to SYMMLQ.
+BiLQ requires support for `adjoint(M)` and `adjoint(N)` if preconditioners are provided.
 
 #### Input arguments
 
@@ -44,6 +46,9 @@ When `A` is Hermitian and `b = c`, BiLQ is equivalent to SYMMLQ.
 
 * `c`: the second initial vector of length `n` required by the Lanczos biorthogonalization process;
 * `transfer_to_bicg`: transfer from the BiLQ point to the BiCG point, when it exists. The transfer is based on the residual norm;
+* `M`: linear operator that models a nonsingular matrix of size `n` used for left preconditioning;
+* `N`: linear operator that models a nonsingular matrix of size `n` used for right preconditioning;
+* `ldiv`: define whether the preconditioners use `ldiv!` or `mul!`;
 * `atol`: absolute stopping tolerance based on the residual norm;
 * `rtol`: relative stopping tolerance based on the residual norm;
 * `itmax`: the maximum number of iterations. If `itmax=0`, the default number of iterations is set to `2n`;
@@ -82,6 +87,9 @@ def_optargs_bilq = (:(x0::AbstractVector),)
 
 def_kwargs_bilq = (:(; c::AbstractVector{FC} = b    ),
                    :(; transfer_to_bicg::Bool = true),
+                   :(; M = I                        ),
+                   :(; N = I                        ),
+                   :(; ldiv::Bool = false           ),
                    :(; atol::T = √eps(T)            ),
                    :(; rtol::T = √eps(T)            ),
                    :(; itmax::Int = 0               ),
@@ -95,7 +103,7 @@ def_kwargs_bilq = mapreduce(extract_parameters, vcat, def_kwargs_bilq)
 
 args_bilq = (:A, :b)
 optargs_bilq = (:x0,)
-kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
+kwargs_bilq = (:c, :transfer_to_bicg, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
 
 @eval begin
   function bilq($(def_args_bilq...), $(def_optargs_bilq...); $(def_kwargs_bilq...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
@@ -131,26 +139,42 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
     length(b) == m || error("Inconsistent problem size")
     (verbose > 0) && @printf(iostream, "BILQ: system of size %d\n", n)
 
+    # Check M = Iₙ and N = Iₙ
+    MisI = (M === I)
+    NisI = (N === I)
+
     # Check type consistency
     eltype(A) == FC || @warn "eltype(A) ≠ $FC. This could lead to errors or additional allocations in operator-vector products."
     ktypeof(b) <: S || error("ktypeof(b) is not a subtype of $S")
     ktypeof(c) <: S || error("ktypeof(c) is not a subtype of $S")
 
-    # Compute the adjoint of A
+    # Compute the adjoint of A, M and N
     Aᴴ = A'
+    Mᴴ = M'
+    Nᴴ = N'
 
     # Set up workspace.
+    allocate_if(!MisI, solver, :t, S, n)
+    allocate_if(!NisI, solver, :s, S, n)
     uₖ₋₁, uₖ, q, vₖ₋₁, vₖ = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ
     p, Δx, x, d̅, stats = solver.p, solver.Δx, solver.x, solver.d̅, solver.stats
     warm_start = solver.warm_start
     rNorms = stats.residuals
     reset!(stats)
     r₀ = warm_start ? q : b
+    Mᴴuₖ = MisI ? uₖ : solver.t
+    t = MisI ? q : solver.t
+    Nvₖ = NisI ? vₖ : solver.s
+    s = NisI ? p : solver.s
 
     if warm_start
       mul!(r₀, A, Δx)
       @kaxpby!(n, one(FC), b, -one(FC), r₀)
     end
+    if !MisI
+      mulorldiv!(solver.t, M, r₀, ldiv)
+      r₀ = solver.t
+    end
 
     # Initial solution x₀ and residual norm ‖r₀‖.
     x .= zero(FC)
@@ -170,10 +194,6 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
     iter = 0
     itmax == 0 && (itmax = 2*n)
 
-    ε = atol + rtol * bNorm
-    (verbose > 0) && @printf(iostream, "%5s  %7s  %5s\n", "k", "‖rₖ‖", "timer")
-    kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %.2fs\n", iter, bNorm, ktimer(start_time))
-
     # Initialize the Lanczos biorthogonalization process.
     cᴴb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
     if cᴴb == 0
@@ -186,6 +206,10 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
       return solver
     end
 
+    ε = atol + rtol * bNorm
+    (verbose > 0) && @printf(iostream, "%5s  %8s  %7s  %5s\n", "k", "αₖ", "‖rₖ‖", "timer")
+    kdisplay(iter, verbose) && @printf(iostream, "%5d  %8.1e  %7.1e  %.2fs\n", iter, cᴴb, bNorm, ktimer(start_time))
+
     βₖ = √(abs(cᴴb))            # β₁γ₁ = cᴴ(b - Ax₀)
     γₖ = cᴴb / βₖ               # β₁γ₁ = cᴴ(b - Ax₀)
     vₖ₋₁ .= zero(FC)            # v₀ = 0
@@ -214,23 +238,30 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
       iter = iter + 1
 
       # Continue the Lanczos biorthogonalization process.
-      # AVₖ  = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
-      # AᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
+      # MANVₖ    = VₖTₖ    + βₖ₊₁vₖ₊₁(eₖ)ᵀ = Vₖ₊₁Tₖ₊₁.ₖ
+      # NᴴAᴴMᴴUₖ = Uₖ(Tₖ)ᴴ + γ̄ₖ₊₁uₖ₊₁(eₖ)ᵀ = Uₖ₊₁(Tₖ.ₖ₊₁)ᴴ
 
-      mul!(q, A , vₖ)  # Forms vₖ₊₁ : q ← Avₖ
-      mul!(p, Aᴴ, uₖ)  # Forms uₖ₊₁ : p ← Aᴴuₖ
+      # Forms vₖ₊₁ : q ← MANvₖ
+      NisI || mulorldiv!(Nvₖ, N, vₖ, ldiv)
+      mul!(t, A, Nvₖ)
+      MisI || mulorldiv!(q, M, t, ldiv)
+
+      # Forms uₖ₊₁ : p ← NᴴAᴴMᴴuₖ
+      MisI || mulorldiv!(Mᴴuₖ, Mᴴ, uₖ, ldiv)
+      mul!(s, Aᴴ, Mᴴuₖ)
+      NisI || mulorldiv!(p, Nᴴ, s, ldiv)
 
       @kaxpy!(n, -γₖ, vₖ₋₁, q)  # q ← q - γₖ * vₖ₋₁
       @kaxpy!(n, -βₖ, uₖ₋₁, p)  # p ← p - β̄ₖ * uₖ₋₁
 
-      αₖ = @kdot(n, uₖ, q)      # αₖ = ⟨uₖ,q⟩
+      αₖ = @kdot(n, uₖ, q)  # αₖ = ⟨uₖ,q⟩
 
-      @kaxpy!(n, -     αₖ , vₖ, q)    # q ← q - αₖ * vₖ
-      @kaxpy!(n, -conj(αₖ), uₖ, p)    # p ← p - ᾱₖ * uₖ
+      @kaxpy!(n, -     αₖ , vₖ, q)  # q ← q - αₖ * vₖ
+      @kaxpy!(n, -conj(αₖ), uₖ, p)  # p ← p - ᾱₖ * uₖ
 
-      pᴴq = @kdot(n, p, q)      # pᴴq  = ⟨p,q⟩
-      βₖ₊₁ = √(abs(pᴴq))        # βₖ₊₁ = √(|pᴴq|)
-      γₖ₊₁ = pᴴq / βₖ₊₁         # γₖ₊₁ = pᴴq / βₖ₊₁
+      pᴴq = @kdot(n, p, q)  # pᴴq  = ⟨p,q⟩
+      βₖ₊₁ = √(abs(pᴴq))    # βₖ₊₁ = √(|pᴴq|)
+      γₖ₊₁ = pᴴq / βₖ₊₁     # γₖ₊₁ = pᴴq / βₖ₊₁
 
       # Update the LQ factorization of Tₖ = L̅ₖQₖ.
       # [ α₁ γ₂ 0  •  •  •  0 ]   [ δ₁   0    •   •   •    •    0   ]
@@ -298,19 +329,19 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
       # Compute d̅ₖ.
       if iter == 1
         # d̅₁ = v₁
-        @. d̅ = vₖ
+        @kcopy!(n, vₖ, d̅)  # d̅ ← vₖ
       else
         # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
         @kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅)
       end
 
       # Compute vₖ₊₁ and uₖ₊₁.
-      @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
-      @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ
+      @kcopy!(n, vₖ, vₖ₋₁)  # vₖ₋₁ ← vₖ
+      @kcopy!(n, uₖ, uₖ₋₁)  # uₖ₋₁ ← uₖ
 
       if pᴴq ≠ 0
-        @. vₖ = q / βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
-        @. uₖ = p / conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
+        vₖ .= q ./ βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
+        uₖ .= p ./ conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
       end
 
       # Compute ⟨vₖ,vₖ₊₁⟩ and ‖vₖ₊₁‖
@@ -353,7 +384,7 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
       breakdown = !solved_lq && !solved_cg && (pᴴq == 0)
       timer = time_ns() - start_time
       overtimed = timer > timemax_ns
-      kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %.2fs\n", iter, rNorm_lq, ktimer(start_time))
+      kdisplay(iter, verbose) && @printf(iostream, "%5d  %8.1e  %7.1e  %.2fs\n", iter, αₖ, rNorm_lq, ktimer(start_time))
     end
     (verbose > 0) && @printf(iostream, "\n")
 
@@ -372,6 +403,10 @@ kwargs_bilq = (:c, :transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose,
     overtimed           && (status = "time limit exceeded")
 
     # Update x
+    if !NisI
+      copyto!(solver.s, x)
+      mulorldiv!(x, N, solver.s, ldiv)
+    end
     warm_start && @kaxpy!(n, one(FC), Δx, x)
     solver.warm_start = false
 

src/bilqr.jl

@@ -316,7 +316,7 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi
         # Compute d̅ₖ.
         if iter == 1
           # d̅₁ = v₁
-          @. d̅ = vₖ
+          @kcopy!(n, vₖ, d̅)  # d̅ ← vₖ
         else
           # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * vₖ
           @kaxpby!(n, -cₖ, vₖ, conj(sₖ), d̅)
@@ -375,22 +375,22 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi
         if iter == 2
           wₖ₋₁ = wₖ₋₂
           @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁)
-          @. wₖ₋₁ = uₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= uₖ₋₁ ./ conj(δₖ₋₁)
         end
         # w₂ = (u₂ - λ̄₁w₁) / δ̄₂
         if iter == 3
           wₖ₋₁ = wₖ₋₃
           @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁)
           @kaxpy!(n, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁)
-          @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁)
         end
         # wₖ₋₁ = (uₖ₋₁ - λ̄ₖ₋₂wₖ₋₂ - ϵ̄ₖ₋₃wₖ₋₃) / δ̄ₖ₋₁
         if iter ≥ 4
           @kscal!(n, -conj(ϵₖ₋₃), wₖ₋₃)
           wₖ₋₁ = wₖ₋₃
           @kaxpy!(n, one(FC), uₖ₋₁, wₖ₋₁)
           @kaxpy!(n, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁)
-          @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁)
         end
 
         if iter ≥ 3
@@ -421,12 +421,12 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi
       end
 
       # Compute vₖ₊₁ and uₖ₊₁.
-      @. vₖ₋₁ = vₖ  # vₖ₋₁ ← vₖ
-      @. uₖ₋₁ = uₖ  # uₖ₋₁ ← uₖ
+      @kcopy!(n, vₖ, vₖ₋₁)  # vₖ₋₁ ← vₖ
+      @kcopy!(n, uₖ, uₖ₋₁)  # uₖ₋₁ ← uₖ
 
       if pᴴq ≠ zero(FC)
-        @. vₖ = q / βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
-        @. uₖ = p / conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
+        vₖ .= q ./ βₖ₊₁        # βₖ₊₁vₖ₊₁ = q
+        uₖ .= p ./ conj(γₖ₊₁)  # γ̄ₖ₊₁uₖ₊₁ = p
       end
 
       # Update ϵₖ₋₃, λₖ₋₂, δbarₖ₋₁, cₖ₋₁, sₖ₋₁, γₖ and βₖ.

src/block_gmres.jl

@@ -148,6 +148,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto
     # Define the blocks D1 and D2
     D1 = view(D, 1:p, :)
     D2 = view(D, p+1:2p, :)
+    trans = FC <: AbstractFloat ? 'T' : 'C'
 
     # Coefficients for mul!
     α = -one(FC)
@@ -263,7 +264,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto
         for i = 1 : inner_iter-1
           D1 .= R[nr+i]
           D2 .= R[nr+i+1]
-          @kormqr!('L', 'T', H[i], τ[i], D)
+          @kormqr!('L', trans, H[i], τ[i], D)
           R[nr+i] .= D1
           R[nr+i+1] .= D2
         end
@@ -276,7 +277,7 @@ kwargs_block_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rto
         # Update Zₖ = (Qₖ)ᴴΓE₁ = (Λ₁, ..., Λₖ, Λbarₖ₊₁)
         D1 .= Z[inner_iter]
         D2 .= zero(FC)
-        @kormqr!('L', 'T', H[inner_iter], τ[inner_iter], D)
+        @kormqr!('L', trans, H[inner_iter], τ[inner_iter], D)
         Z[inner_iter] .= D1
 
         # Update residual norm estimate.

src/block_krylov_utils.jl

@@ -164,7 +164,13 @@ function copy_triangle(Q::AbstractMatrix{FC}, R::AbstractMatrix{FC}, k::Int) whe
       end
     end
   else
-    copytrito!(R, Q, 'U')
+    mR, nR = size(R)
+    mQ, nQ = size(Q)
+    if (mR == mQ) && (nR == nQ)
+      copytrito!(R, Q, 'U')
+    else
+      copytrito!(R, view(Q, 1:k, 1:k), 'U')
+    end
   end
   return R
 end

src/cg_lanczos.jl

@@ -218,9 +218,9 @@ kwargs_cg_lanczos = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :timemax
       @kaxpy!(n, -δ, Mv, Mv_next)        # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ
       if iter > 0
         @kaxpy!(n, -β, Mv_prev, Mv_next) # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁
-        @. Mv_prev = Mv                  # Mvₖ₋₁ ← Mvₖ
+        @kcopy!(n, Mv, Mv_prev)          # Mvₖ₋₁ ← Mvₖ
       end
-      @. Mv = Mv_next                      # Mvₖ ← Mvₖ₊₁
+      @kcopy!(n, Mv_next, Mv)              # Mvₖ ← Mvₖ₊₁
       MisI || mulorldiv!(v, M, Mv, ldiv)   # vₖ₊₁ = M⁻¹ * Mvₖ₊₁
       β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁
       @kscal!(n, one(FC) / β, v)           # vₖ₊₁  ←  vₖ₊₁ / βₖ₊₁

src/cg_lanczos_shift.jl

@@ -211,9 +211,9 @@ kwargs_cg_lanczos_shift = (:M, :ldiv, :check_curvature, :atol, :rtol, :itmax, :t
       @kaxpy!(n, -δ, Mv, Mv_next)          # Mvₖ₊₁ ← Mvₖ₊₁ - δₖMvₖ
       if iter > 0
         @kaxpy!(n, -β, Mv_prev, Mv_next)   # Mvₖ₊₁ ← Mvₖ₊₁ - βₖMvₖ₋₁
-        @. Mv_prev = Mv                    # Mvₖ₋₁ ← Mvₖ
+        @kcopy!(n, Mv, Mv_prev)            # Mvₖ₋₁ ← Mvₖ
       end
-      @. Mv = Mv_next                      # Mvₖ ← Mvₖ₊₁
+      @kcopy!(n, Mv_next, Mv)              # Mvₖ ← Mvₖ₊₁
       MisI || mulorldiv!(v, M, Mv, ldiv)   # vₖ₊₁ = M⁻¹ * Mvₖ₊₁
       β = sqrt(@kdotr(n, v, Mv))           # βₖ₊₁ = vₖ₊₁ᴴ M vₖ₊₁
       @kscal!(n, one(FC) / β, v)           # vₖ₊₁  ←  vₖ₊₁ / βₖ₊₁

src/fgmres.jl

@@ -241,7 +241,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i
       # Initial ζ₁ and V₁
       β = @knrm2(n, r₀)
       z[1] = β
-      @. V[1] = r₀ / rNorm
+      V[1] .= r₀ ./ rNorm
 
       npass = npass + 1
       solver.inner_iter = 0
@@ -332,7 +332,7 @@ kwargs_fgmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :i
             push!(V, S(undef, n))
             push!(z, zero(FC))
           end
-          @. V[inner_iter+1] = q / Hbis  # hₖ₊₁.ₖvₖ₊₁ = q
+          V[inner_iter+1] .= q ./ Hbis  # hₖ₊₁.ₖvₖ₊₁ = q
           z[inner_iter+1] = ζₖ₊₁
         end
       end

src/fom.jl

@@ -233,7 +233,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma
       # Initial ζ₁ and V₁
       β = @knrm2(n, r₀)
       z[1] = β
-      @. V[1] = r₀ / rNorm
+      V[1] .= r₀ ./ rNorm
 
       npass = npass + 1
       inner_iter = 0
@@ -314,7 +314,7 @@ kwargs_fom = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :itma
           if !restart && (inner_iter ≥ mem)
             push!(V, S(undef, n))
           end
-          @. V[inner_iter+1] = q / Hbis  # hₖ₊₁.ₖvₖ₊₁ = q
+          V[inner_iter+1] .= q ./ Hbis  # hₖ₊₁.ₖvₖ₊₁ = q
         end
       end