diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 120f1f736..358eb9cac 100644
--- a/.buildkite/pipeline.yml
+++ b/.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
diff --git a/docs/src/preconditioners.md b/docs/src/preconditioners.md
index 22e24ce90..38da24cc5 100644
--- a/docs/src/preconditioners.md
+++ b/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.
 
diff --git a/src/bilq.jl b/src/bilq.jl
index 2e8823e93..5bda9802a 100644
--- a/src/bilq.jl
+++ b/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
 
diff --git a/src/bilqr.jl b/src/bilqr.jl
index 486ccceec..1abc60bce 100644
--- a/src/bilqr.jl
+++ b/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,14 +375,14 @@ 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
@@ -390,7 +390,7 @@ kwargs_bilqr = (:transfer_to_bicg, :atol, :rtol, :itmax, :timemax, :verbose, :hi
           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 βₖ.
diff --git a/src/block_gmres.jl b/src/block_gmres.jl
index c1e04164d..ea0b4b989 100644
--- a/src/block_gmres.jl
+++ b/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.
diff --git a/src/block_krylov_utils.jl b/src/block_krylov_utils.jl
index 68ba46458..9f1a9616d 100644
--- a/src/block_krylov_utils.jl
+++ b/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
diff --git a/src/cg_lanczos.jl b/src/cg_lanczos.jl
index 2c5d72a64..778a73d6f 100644
--- a/src/cg_lanczos.jl
+++ b/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ₖ₊₁ / βₖ₊₁
diff --git a/src/cg_lanczos_shift.jl b/src/cg_lanczos_shift.jl
index b523e5cc3..e47c882e6 100644
--- a/src/cg_lanczos_shift.jl
+++ b/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ₖ₊₁ / βₖ₊₁
diff --git a/src/fgmres.jl b/src/fgmres.jl
index 1a68aac6c..0972f9274 100644
--- a/src/fgmres.jl
+++ b/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
diff --git a/src/fom.jl b/src/fom.jl
index 351fb246f..c27485f2f 100644
--- a/src/fom.jl
+++ b/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
 
diff --git a/src/gmres.jl b/src/gmres.jl
index 7ee6e2341..a172be6d9 100644
--- a/src/gmres.jl
+++ b/src/gmres.jl
@@ -235,7 +235,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it
       # Initial ζ₁ and V₁
       β = @knrm2(n, r₀)
       z[1] = β
-      @. V[1] = r₀ / rNorm
+      V[1] .= r₀ ./ rNorm
 
       npass = npass + 1
       solver.inner_iter = 0
@@ -325,7 +325,7 @@ kwargs_gmres = (:M, :N, :ldiv, :restart, :reorthogonalization, :atol, :rtol, :it
             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
diff --git a/src/gpmr.jl b/src/gpmr.jl
index 1049c3b50..b47d9d32e 100644
--- a/src/gpmr.jl
+++ b/src/gpmr.jl
@@ -261,12 +261,12 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato
     # βv₁ = Cb
     β = @knrm2(m, b₀)
     β ≠ 0 || error("b must be nonzero")
-    @. V[1] = b₀ / β
+    V[1] .= b₀ ./ β
 
     # γu₁ = Dc
     γ = @knrm2(n, c₀)
     γ ≠ 0 || error("c must be nonzero")
-    @. U[1] = c₀ / γ
+    U[1] .= c₀ ./ γ
 
     # Compute ‖r₀‖² = γ² + β²
     rNorm = sqrt(γ^2 + β^2)
@@ -479,7 +479,7 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato
 
         # hₖ₊₁.ₖ ≠ 0
         if Haux > btol
-          @. V[k+1] = q / Haux  # hₖ₊₁.ₖvₖ₊₁ = q
+          V[k+1] .= q ./ Haux  # hₖ₊₁.ₖvₖ₊₁ = q
         else
           # Breakdown -- hₖ₊₁.ₖ = ‖q‖₂ = 0 and Auₖ ∈ Span{v₁, ..., vₖ}
           V[k+1] .= zero(FC)  # vₖ₊₁ = 0 such that vₖ₊₁ ⊥ Span{v₁, ..., vₖ}
@@ -487,7 +487,7 @@ kwargs_gpmr = (:C, :D, :E, :F, :ldiv, :gsp, :λ, :μ, :reorthogonalization, :ato
 
         # fₖ₊₁.ₖ ≠ 0
         if Faux > btol
-          @. U[k+1] = p / Faux  # fₖ₊₁.ₖuₖ₊₁ = p
+          U[k+1] .= p ./ Faux  # fₖ₊₁.ₖuₖ₊₁ = p
         else
           # Breakdown -- fₖ₊₁.ₖ = ‖p‖₂ = 0 and Bvₖ ∈ Span{u₁, ..., uₖ}
           U[k+1] .= zero(FC)  # uₖ₊₁ = 0 such that uₖ₊₁ ⊥ Span{u₁, ..., uₖ}
diff --git a/src/krylov_solvers.jl b/src/krylov_solvers.jl
index 13c98f823..d9f9fd9b1 100644
--- a/src/krylov_solvers.jl
+++ b/src/krylov_solvers.jl
@@ -1002,6 +1002,8 @@ mutable struct BilqSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   Δx         :: S
   x          :: S
   d̅          :: S
+  t          :: S
+  s          :: S
   warm_start :: Bool
   stats      :: SimpleStats{T}
 end
@@ -1018,8 +1020,10 @@ function BilqSolver(m, n, S)
   Δx   = S(undef, 0)
   x    = S(undef, n)
   d̅    = S(undef, n)
+  t    = S(undef, 0)
+  s    = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown")
-  solver = BilqSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, d̅, false, stats)
+  solver = BilqSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, d̅, t, s, false, stats)
   return solver
 end
 
@@ -1052,6 +1056,8 @@ mutable struct QmrSolver{T,FC,S} <: KrylovSolver{T,FC,S}
   x          :: S
   wₖ₋₂       :: S
   wₖ₋₁       :: S
+  t          :: S
+  s          :: S
   warm_start :: Bool
   stats      :: SimpleStats{T}
 end
@@ -1069,8 +1075,10 @@ function QmrSolver(m, n, S)
   x    = S(undef, n)
   wₖ₋₂ = S(undef, n)
   wₖ₋₁ = S(undef, n)
+  t    = S(undef, 0)
+  s    = S(undef, 0)
   stats = SimpleStats(0, false, false, T[], T[], T[], 0.0, "unknown")
-  solver = QmrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, wₖ₋₂, wₖ₋₁, false, stats)
+  solver = QmrSolver{T,FC,S}(m, n, uₖ₋₁, uₖ, q, vₖ₋₁, vₖ, p, Δx, x, wₖ₋₂, wₖ₋₁, t, s, false, stats)
   return solver
 end
 
diff --git a/src/minres.jl b/src/minres.jl
index d6ca85d8c..2e0dd114b 100644
--- a/src/minres.jl
+++ b/src/minres.jl
@@ -276,8 +276,8 @@ kwargs_minres = (:M, :ldiv, :λ, :atol, :rtol, :etol, :conlim, :itmax, :timemax,
       end
       @kaxpy!(n, one(FC) / β, v, w)
 
-      @. r1 = r2
-      @. r2 = y
+      @kcopy!(n, r2, r1)  # r1 ← r2
+      @kcopy!(n, y, r2)   # r2 ← y
       MisI || mulorldiv!(v, M, r2, ldiv)
       oldβ = β
       β = @kdotr(n, r2, v)
diff --git a/src/minres_qlp.jl b/src/minres_qlp.jl
index 5bc3399eb..9cb569f6a 100644
--- a/src/minres_qlp.jl
+++ b/src/minres_qlp.jl
@@ -369,12 +369,12 @@ kwargs_minres_qlp = (:M, :ldiv, :λ, :atol, :rtol, :Artol, :itmax, :timemax, :ve
       # Compute directions wₖ₋₂, ẘₖ₋₁ and w̄ₖ, last columns of Wₖ = Vₖ(Pₖ)ᴴ
       if iter == 1
         # w̅₁ = v₁
-        @. wₖ = vₖ
+        @kcopy!(n, vₖ, wₖ)
       elseif iter == 2
         # [w̅ₖ₋₁ vₖ] [cpₖ  spₖ] = [ẘₖ₋₁ w̅ₖ] ⟷ ẘₖ₋₁ = cpₖ * w̅ₖ₋₁ + spₖ * vₖ
         #           [spₖ -cpₖ]             ⟷ w̅ₖ   = spₖ * w̅ₖ₋₁ - cpₖ * vₖ
         @kswap(wₖ₋₁, wₖ)
-        @. wₖ = spₖ * wₖ₋₁ - cpₖ * vₖ
+        wₖ .= spₖ .* wₖ₋₁ .- cpₖ .* vₖ
         @kaxpby!(n, spₖ, vₖ, cpₖ, wₖ₋₁)
       else
         # [ẘₖ₋₂ w̄ₖ₋₁ vₖ] [cpₖ  0   spₖ] [1   0    0 ] = [wₖ₋₂ ẘₖ₋₁ w̄ₖ] ⟷ wₖ₋₂ = cpₖ * ẘₖ₋₂ + spₖ * vₖ
diff --git a/src/qmr.jl b/src/qmr.jl
index 995392f0c..2a1ad578c 100644
--- a/src/qmr.jl
+++ b/src/qmr.jl
@@ -22,7 +22,7 @@ export qmr, qmr!
 
 """
     (x, stats) = qmr(A, b::AbstractVector{FC};
-                     c::AbstractVector{FC}=b, atol::T=√eps(T),
+                     c::AbstractVector{FC}=b, 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)
 
@@ -38,6 +38,7 @@ Solve the square linear system Ax = b of size n using QMR.
 QMR 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`, QMR is equivalent to MINRES.
+QMR requires support for `adjoint(M)` and `adjoint(N)` if preconditioners are provided.
 
 #### Input arguments
 
@@ -51,6 +52,9 @@ When `A` is Hermitian and `b = c`, QMR is equivalent to MINRES.
 #### Keyword arguments
 
 * `c`: the second initial vector of length `n` required by the Lanczos biorthogonalization process;
+* `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`;
@@ -89,6 +93,9 @@ def_args_qmr = (:(A                    ),
 def_optargs_qmr = (:(x0::AbstractVector),)
 
 def_kwargs_qmr = (:(; c::AbstractVector{FC} = b ),
+                  :(; M = I                     ),
+                  :(; N = I                     ),
+                  :(; ldiv::Bool = false        ),
                   :(; atol::T = √eps(T)         ),
                   :(; rtol::T = √eps(T)         ),
                   :(; itmax::Int = 0            ),
@@ -102,7 +109,7 @@ def_kwargs_qmr = mapreduce(extract_parameters, vcat, def_kwargs_qmr)
 
 args_qmr = (:A, :b)
 optargs_qmr = (:x0,)
-kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
+kwargs_qmr = (:c, :M, :N, :ldiv, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback, :iostream)
 
 @eval begin
   function qmr($(def_args_qmr...), $(def_optargs_qmr...); $(def_kwargs_qmr...)) where {T <: AbstractFloat, FC <: FloatOrComplex{T}}
@@ -138,26 +145,42 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
     length(b) == m || error("Inconsistent problem size")
     (verbose > 0) && @printf(iostream, "QMR: 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ₖ, p = solver.uₖ₋₁, solver.uₖ, solver.q, solver.vₖ₋₁, solver.vₖ, solver.p
     Δx, x, wₖ₋₂, wₖ₋₁, stats = solver.Δx, solver.x, solver.wₖ₋₂, solver.wₖ₋₁, 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)
@@ -177,10 +200,6 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
     iter = 0
     itmax == 0 && (itmax = 2*n)
 
-    ε = atol + rtol * rNorm
-    (verbose > 0) && @printf(iostream, "%5s  %7s  %5s\n", "k", "‖rₖ‖", "timer")
-    kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %.2fs\n", iter, rNorm, ktimer(start_time))
-
     # Initialize the Lanczos biorthogonalization process.
     cᴴb = @kdot(n, c, r₀)  # ⟨c,r₀⟩
     if cᴴb == 0
@@ -193,6 +212,10 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       return solver
     end
 
+    ε = atol + rtol * rNorm
+    (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, rNorm, ktimer(start_time))
+
     βₖ = √(abs(cᴴb))             # β₁γ₁ = cᴴ(b - Ax₀)
     γₖ = cᴴb / βₖ                # β₁γ₁ = cᴴ(b - Ax₀)
     vₖ₋₁ .= zero(FC)             # v₀ = 0
@@ -219,23 +242,30 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       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 QR factorization of Tₖ₊₁.ₖ = Qₖ [ Rₖ ].
       #                                            [ Oᵀ ]
@@ -296,14 +326,14 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       if iter == 1
         wₖ = wₖ₋₁
         @kaxpy!(n, one(FC), vₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
       # w₂ = (v₂ - λ₁w₁) / δ₂
       if iter == 2
         wₖ = wₖ₋₂
         @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ)
         @kaxpy!(n, one(FC), vₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
       # wₖ = (vₖ - λₖ₋₁wₖ₋₁ - ϵₖ₋₂wₖ₋₂) / δₖ
       if iter ≥ 3
@@ -311,7 +341,7 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
         wₖ = wₖ₋₂
         @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ)
         @kaxpy!(n, one(FC), vₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
 
       # Compute solution xₖ.
@@ -319,12 +349,12 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       @kaxpy!(n, ζₖ, wₖ, x)
 
       # 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
 
       # Compute τₖ₊₁ = τₖ + ‖vₖ₊₁‖²
@@ -357,7 +387,7 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       breakdown = !solved && (pᴴq == 0)
       timer = time_ns() - start_time
       overtimed = timer > timemax_ns
-      kdisplay(iter, verbose) && @printf(iostream, "%5d  %7.1e  %.2fs\n", iter, rNorm, ktimer(start_time))
+      kdisplay(iter, verbose) && @printf(iostream, "%5d  %8.1e  %7.1e  %.2fs\n", iter, αₖ, rNorm, ktimer(start_time))
     end
     (verbose > 0) && @printf(iostream, "\n")
 
@@ -369,6 +399,10 @@ kwargs_qmr = (:c, :atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
     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
 
diff --git a/src/symmlq.jl b/src/symmlq.jl
index 604698525..bfc866944 100644
--- a/src/symmlq.jl
+++ b/src/symmlq.jl
@@ -303,9 +303,9 @@ kwargs_symmlq = (:M, :ldiv, :transfer_to_cg, :λ, :λest, :atol, :rtol, :etol, :
       mul!(Mv_next, A, v)
       α = @kdotr(m, v, Mv_next) + λ
       @kaxpy!(m, -oldβ, Mvold, Mv_next)
-      @. Mvold = Mv
+      @kcopy!(m, Mv, Mvold)  # Mvold ← Mv
       @kaxpy!(m, -α, Mv, Mv_next)
-      @. Mv = Mv_next
+      @kcopy!(m, Mv_next, Mv)  # Mv ← Mv_next
       MisI || mulorldiv!(v, M, Mv, ldiv)
       β = @kdotr(m, v, Mv)
       β < 0 && error("Preconditioner is not positive definite")
diff --git a/src/tricg.jl b/src/tricg.jl
index 8250e6dfc..99d714bdd 100644
--- a/src/tricg.jl
+++ b/src/tricg.jl
@@ -369,22 +369,22 @@ kwargs_tricg = (:M, :N, :ldiv, :spd, :snd, :flip, :τ, :ν, :atol, :rtol, :itmax
       if iter == 1
         # [ 1  0 ] [ gx₁ gy₁ ] = [ v₁ 0  ]
         # [ δ̄₁ 1 ] [ gx₂ gy₂ ]   [ 0  u₁ ]
-        @. gx₂ₖ₋₁ = vₖ
-        @. gx₂ₖ   = - conj(δₖ) * gx₂ₖ₋₁
-        @. gy₂ₖ   = uₖ
+        @kcopy!(m, vₖ, gx₂ₖ₋₁)  # gx₂ₖ₋₁ ← vₖ
+        gx₂ₖ .= -conj(δₖ) .* gx₂ₖ₋₁
+        @kcopy!(n, uₖ, gy₂ₖ)  # gy₂ₖ ← uₖ
       else
         # [ 0  σ̄ₖ 1  0 ] [ gx₂ₖ₋₃ gy₂ₖ₋₃ ] = [ vₖ 0  ]
         # [ η̄ₖ λ̄ₖ δ̄ₖ 1 ] [ gx₂ₖ₋₂ gy₂ₖ₋₂ ]   [ 0  uₖ ]
         #                [ gx₂ₖ₋₁ gy₂ₖ₋₁ ]
         #                [ gx₂ₖ   gy₂ₖ   ]
-        @. gx₂ₖ₋₁ = conj(ηₖ) * gx₂ₖ₋₁ + conj(λₖ) * gx₂ₖ
-        @. gy₂ₖ₋₁ = conj(ηₖ) * gy₂ₖ₋₁ + conj(λₖ) * gy₂ₖ
+        gx₂ₖ₋₁ .= conj(ηₖ) .* gx₂ₖ₋₁ .+ conj(λₖ) .* gx₂ₖ
+        gy₂ₖ₋₁ .= conj(ηₖ) .* gy₂ₖ₋₁ .+ conj(λₖ) .* gy₂ₖ
 
-        @. gx₂ₖ = vₖ - conj(σₖ) * gx₂ₖ
-        @. gy₂ₖ =    - conj(σₖ) * gy₂ₖ
+        gx₂ₖ .= vₖ .- conj(σₖ) .* gx₂ₖ
+        gy₂ₖ .=    .- conj(σₖ) .* gy₂ₖ
 
-        @. gx₂ₖ₋₁ =    - gx₂ₖ₋₁ - conj(δₖ) * gx₂ₖ
-        @. gy₂ₖ₋₁ = uₖ - gy₂ₖ₋₁ - conj(δₖ) * gy₂ₖ
+        gx₂ₖ₋₁ .=    .- gx₂ₖ₋₁ .- conj(δₖ) .* gx₂ₖ
+        gy₂ₖ₋₁ .= uₖ .- gy₂ₖ₋₁ .- conj(δₖ) .* gy₂ₖ
 
         # g₂ₖ₋₃ == g₂ₖ and g₂ₖ₋₂ == g₂ₖ₋₁
         @kswap(gx₂ₖ₋₁, gx₂ₖ)
diff --git a/src/trilqr.jl b/src/trilqr.jl
index 2b584c216..042db2604 100644
--- a/src/trilqr.jl
+++ b/src/trilqr.jl
@@ -300,7 +300,7 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
         # Compute d̅ₖ.
         if iter == 1
           # d̅₁ = u₁
-          @. d̅ = uₖ
+          @kcopy!(n, uₖ, d̅)  # d̅ ← uₖ
         else
           # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ
           @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅)
@@ -351,14 +351,14 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
         if iter == 2
           wₖ₋₁ = wₖ₋₂
           @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁)
-          @. wₖ₋₁ = vₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= vₖ₋₁ ./ conj(δₖ₋₁)
         end
         # w₂ = (v₂ - λ̄₁w₁) / δ̄₂
         if iter == 3
           wₖ₋₁ = wₖ₋₃
           @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁)
           @kaxpy!(m, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁)
-          @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁)
         end
         # wₖ₋₁ = (vₖ₋₁ - λ̄ₖ₋₂wₖ₋₂ - ϵ̄ₖ₋₃wₖ₋₃) / δ̄ₖ₋₁
         if iter ≥ 4
@@ -366,7 +366,7 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
           wₖ₋₁ = wₖ₋₃
           @kaxpy!(m, one(FC), vₖ₋₁, wₖ₋₁)
           @kaxpy!(m, -conj(λₖ₋₂), wₖ₋₂, wₖ₋₁)
-          @. wₖ₋₁ = wₖ₋₁ / conj(δₖ₋₁)
+          wₖ₋₁ .= wₖ₋₁ ./ conj(δₖ₋₁)
         end
 
         if iter ≥ 3
@@ -399,14 +399,14 @@ kwargs_trilqr = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
       end
 
       # Compute uₖ₊₁ and uₖ₊₁.
-      @. vₖ₋₁ = vₖ  # vₖ₋₁ ← vₖ
-      @. uₖ₋₁ = uₖ  # uₖ₋₁ ← uₖ
+      @kcopy!(m, vₖ, vₖ₋₁)  # vₖ₋₁ ← vₖ
+      @kcopy!(n, uₖ, uₖ₋₁)  # uₖ₋₁ ← uₖ
 
       if βₖ₊₁ ≠ zero(T)
-        @. vₖ = q / βₖ₊₁  # βₖ₊₁vₖ₊₁ = q
+        vₖ .= q ./ βₖ₊₁  # βₖ₊₁vₖ₊₁ = q
       end
       if γₖ₊₁ ≠ zero(T)
-        @. uₖ = p / γₖ₊₁  # γₖ₊₁uₖ₊₁ = p
+        uₖ .= p ./ γₖ₊₁  # γₖ₊₁uₖ₊₁ = p
       end
 
       # Update ϵₖ₋₃, λₖ₋₂, δbarₖ₋₁, cₖ₋₁, sₖ₋₁, γₖ and βₖ.
diff --git a/src/trimr.jl b/src/trimr.jl
index ae61b785a..780265e4a 100644
--- a/src/trimr.jl
+++ b/src/trimr.jl
@@ -447,9 +447,9 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
       if iter == 1
         # [ δ₁  0  ] [ gx₁ gy₁ ] = [ v₁ 0  ]
         # [ σ₁  δ₂ ] [ gx₂ gy₂ ]   [ 0  u₁ ]
-        @. gx₂ₖ₋₁ = vₖ / δ₂ₖ₋₁
-        @. gx₂ₖ   = - σ₂ₖ₋₁ / δ₂ₖ * gx₂ₖ₋₁
-        @. gy₂ₖ   = uₖ / δ₂ₖ
+        gx₂ₖ₋₁ .= vₖ ./ δ₂ₖ₋₁
+        gx₂ₖ .= -(σ₂ₖ₋₁ / δ₂ₖ) .* gx₂ₖ₋₁
+        gy₂ₖ .= uₖ ./ δ₂ₖ
       elseif iter == 2
         # [ η₁ σ₂ δ₃ 0  ] [ gx₁ gy₁ ] = [ v₂ 0  ]
         # [ λ₁ η₂ σ₃ δ₄ ] [ gx₂ gy₂ ]   [ 0  u₂ ]
@@ -458,23 +458,23 @@ kwargs_trimr = (:M, :N, :ldiv, :spd, :snd, :flip, :sp, :τ, :ν, :atol, :rtol, :
         @kswap(gx₂ₖ₋₃, gx₂ₖ₋₁)
         @kswap(gx₂ₖ₋₂, gx₂ₖ)
         @kswap(gy₂ₖ₋₂, gy₂ₖ)
-        @. gx₂ₖ₋₁ = (vₖ - η₂ₖ₋₃ * gx₂ₖ₋₃ - σ₂ₖ₋₂ * gx₂ₖ₋₂                 ) / δ₂ₖ₋₁
-        @. gx₂ₖ   = (   - λ₂ₖ₋₃ * gx₂ₖ₋₃ - η₂ₖ₋₂ * gx₂ₖ₋₂ - σ₂ₖ₋₁ * gx₂ₖ₋₁) / δ₂ₖ
-        @. gy₂ₖ₋₁ = (   - η₂ₖ₋₃ * gy₂ₖ₋₃ - σ₂ₖ₋₂ * gy₂ₖ₋₂                 ) / δ₂ₖ₋₁
-        @. gy₂ₖ   = (uₖ - λ₂ₖ₋₃ * gy₂ₖ₋₃ - η₂ₖ₋₂ * gy₂ₖ₋₂ - σ₂ₖ₋₁ * gy₂ₖ₋₁) / δ₂ₖ
+        gx₂ₖ₋₁ .= (vₖ .- η₂ₖ₋₃ .* gx₂ₖ₋₃ .- σ₂ₖ₋₂ .* gx₂ₖ₋₂                   ) ./ δ₂ₖ₋₁
+        gx₂ₖ   .= (   .- λ₂ₖ₋₃ .* gx₂ₖ₋₃ .- η₂ₖ₋₂ .* gx₂ₖ₋₂ .- σ₂ₖ₋₁ .* gx₂ₖ₋₁) ./ δ₂ₖ
+        gy₂ₖ₋₁ .= (   .- η₂ₖ₋₃ .* gy₂ₖ₋₃ .- σ₂ₖ₋₂ .* gy₂ₖ₋₂                   ) ./ δ₂ₖ₋₁
+        gy₂ₖ   .= (uₖ .- λ₂ₖ₋₃ .* gy₂ₖ₋₃ .- η₂ₖ₋₂ .* gy₂ₖ₋₂ .- σ₂ₖ₋₁ .* gy₂ₖ₋₁) ./ δ₂ₖ
       else
         # μ₂ₖ₋₅ * gx₂ₖ₋₅ + λ₂ₖ₋₄ * gx₂ₖ₋₄ + η₂ₖ₋₃ * gx₂ₖ₋₃ + σ₂ₖ₋₂ * gx₂ₖ₋₂ + δ₂ₖ₋₁ * gx₂ₖ₋₁              = vₖ
         #                  μ₂ₖ₋₄ * gx₂ₖ₋₄ + λ₂ₖ₋₃ * gx₂ₖ₋₃ + η₂ₖ₋₂ * gx₂ₖ₋₂ + σ₂ₖ₋₁ * gx₂ₖ₋₁ + δ₂ₖ * gx₂ₖ = 0
         g₂ₖ₋₁ = g₂ₖ₋₅ = gx₂ₖ₋₃; g₂ₖ = g₂ₖ₋₄ = gx₂ₖ₋₂; g₂ₖ₋₃ = gx₂ₖ₋₁; g₂ₖ₋₂ = gx₂ₖ
-        @. g₂ₖ₋₁ = (vₖ - μ₂ₖ₋₅ * g₂ₖ₋₅ - λ₂ₖ₋₄ * g₂ₖ₋₄ - η₂ₖ₋₃ * g₂ₖ₋₃ - σ₂ₖ₋₂ * g₂ₖ₋₂                ) / δ₂ₖ₋₁
-        @. g₂ₖ   = (                   - μ₂ₖ₋₄ * g₂ₖ₋₄ - λ₂ₖ₋₃ * g₂ₖ₋₃ - η₂ₖ₋₂ * g₂ₖ₋₂ - σ₂ₖ₋₁ * g₂ₖ₋₁) / δ₂ₖ
+        g₂ₖ₋₁ .= (vₖ .- μ₂ₖ₋₅ .* g₂ₖ₋₅ .- λ₂ₖ₋₄ .* g₂ₖ₋₄ .- η₂ₖ₋₃ .* g₂ₖ₋₃ .- σ₂ₖ₋₂ .* g₂ₖ₋₂                  ) ./ δ₂ₖ₋₁
+        g₂ₖ   .= (                     .- μ₂ₖ₋₄ .* g₂ₖ₋₄ .- λ₂ₖ₋₃ .* g₂ₖ₋₃ .- η₂ₖ₋₂ .* g₂ₖ₋₂ .- σ₂ₖ₋₁ .* g₂ₖ₋₁) ./ δ₂ₖ
         @kswap(gx₂ₖ₋₃, gx₂ₖ₋₁)
         @kswap(gx₂ₖ₋₂, gx₂ₖ)
         # μ₂ₖ₋₅ * gy₂ₖ₋₅ + λ₂ₖ₋₄ * gy₂ₖ₋₄ + η₂ₖ₋₃ * gy₂ₖ₋₃ + σ₂ₖ₋₂ * gy₂ₖ₋₂ + δ₂ₖ₋₁ * gy₂ₖ₋₁              = 0
         #                  μ₂ₖ₋₄ * gy₂ₖ₋₄ + λ₂ₖ₋₃ * gy₂ₖ₋₃ + η₂ₖ₋₂ * gy₂ₖ₋₂ + σ₂ₖ₋₁ * gy₂ₖ₋₁ + δ₂ₖ * gy₂ₖ = uₖ
         g₂ₖ₋₁ = g₂ₖ₋₅ = gy₂ₖ₋₃; g₂ₖ = g₂ₖ₋₄ = gy₂ₖ₋₂; g₂ₖ₋₃ = gy₂ₖ₋₁; g₂ₖ₋₂ = gy₂ₖ
-        @. g₂ₖ₋₁ = (     - μ₂ₖ₋₅ * g₂ₖ₋₅ - λ₂ₖ₋₄ * g₂ₖ₋₄ - η₂ₖ₋₃ * g₂ₖ₋₃ - σ₂ₖ₋₂ * g₂ₖ₋₂                ) / δ₂ₖ₋₁
-        @. g₂ₖ   = (uₖ                   - μ₂ₖ₋₄ * g₂ₖ₋₄ - λ₂ₖ₋₃ * g₂ₖ₋₃ - η₂ₖ₋₂ * g₂ₖ₋₂ - σ₂ₖ₋₁ * g₂ₖ₋₁) / δ₂ₖ
+        g₂ₖ₋₁ .= (     .- μ₂ₖ₋₅ .* g₂ₖ₋₅ .- λ₂ₖ₋₄ .* g₂ₖ₋₄ .- η₂ₖ₋₃ .* g₂ₖ₋₃ .- σ₂ₖ₋₂ .* g₂ₖ₋₂                  ) ./ δ₂ₖ₋₁
+        g₂ₖ   .= (uₖ                     .- μ₂ₖ₋₄ .* g₂ₖ₋₄ .- λ₂ₖ₋₃ .* g₂ₖ₋₃ .- η₂ₖ₋₂ .* g₂ₖ₋₂ .- σ₂ₖ₋₁ .* g₂ₖ₋₁) ./ δ₂ₖ
         @kswap(gy₂ₖ₋₃, gy₂ₖ₋₁)
         @kswap(gy₂ₖ₋₂, gy₂ₖ)
       end
diff --git a/src/usymlq.jl b/src/usymlq.jl
index b80f0a622..73d6c7114 100644
--- a/src/usymlq.jl
+++ b/src/usymlq.jl
@@ -38,7 +38,8 @@ USYMLQ determines the least-norm solution of the consistent linear system Ax = b
 USYMLQ is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
 The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`.
 The error norm ‖x - x*‖ monotonously decreases in USYMLQ.
-It's considered as a generalization of SYMMLQ.
+When `A` is Hermitian and `b = c`, USYMLQ is equivalent to SYMMLQ.
+USYMLQ is considered as a generalization of SYMMLQ.
 
 It can also be applied to under-determined and over-determined problems.
 In all cases, problems must be consistent.
@@ -295,21 +296,21 @@ kwargs_usymlq = (:transfer_to_usymcg, :atol, :rtol, :itmax, :timemax, :verbose,
       # Compute d̅ₖ.
       if iter == 1
         # d̅₁ = u₁
-        @. d̅ = uₖ
+        @kcopy!(n, uₖ, d̅)  # d̅ ← vₖ
       else
         # d̅ₖ = s̄ₖ * d̅ₖ₋₁ - cₖ * uₖ
         @kaxpby!(n, -cₖ, uₖ, conj(sₖ), d̅)
       end
 
       # Compute uₖ₊₁ and uₖ₊₁.
-      @. vₖ₋₁ = vₖ  # vₖ₋₁ ← vₖ
-      @. uₖ₋₁ = uₖ  # uₖ₋₁ ← uₖ
+      @kcopy!(m, vₖ, vₖ₋₁)  # vₖ₋₁ ← vₖ
+      @kcopy!(n, uₖ, uₖ₋₁)  # uₖ₋₁ ← uₖ
 
       if βₖ₊₁ ≠ zero(T)
-        @. vₖ = q / βₖ₊₁  # βₖ₊₁vₖ₊₁ = q
+        vₖ .= q ./ βₖ₊₁  # βₖ₊₁vₖ₊₁ = q
       end
       if γₖ₊₁ ≠ zero(T)
-        @. uₖ = p / γₖ₊₁  # γₖ₊₁uₖ₊₁ = p
+        uₖ .= p ./ γₖ₊₁  # γₖ₊₁uₖ₊₁ = p
       end
 
       # Compute USYMLQ residual norm
diff --git a/src/usymqr.jl b/src/usymqr.jl
index 0aae23335..1d7488e36 100644
--- a/src/usymqr.jl
+++ b/src/usymqr.jl
@@ -38,7 +38,8 @@ USYMQR solves Ax = b if it is consistent.
 USYMQR is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors `b` and `c`.
 The vector `c` is only used to initialize the process and a default value can be `b` or `Aᴴb` depending on the shape of `A`.
 The residual norm ‖b - Ax‖ monotonously decreases in USYMQR.
-It's considered as a generalization of MINRES.
+When `A` is Hermitian and `b = c`, QMR is equivalent to MINRES.
+USYMQR is considered as a generalization of MINRES.
 
 It can also be applied to under-determined and over-determined problems.
 USYMQR finds the minimum-norm solution if problems are inconsistent.
@@ -278,14 +279,14 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       if iter == 1
         wₖ = wₖ₋₁
         @kaxpy!(n, one(FC), uₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
       # w₂ = (u₂ - λ₁w₁) / δ₂
       if iter == 2
         wₖ = wₖ₋₂
         @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ)
         @kaxpy!(n, one(FC), uₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
       # wₖ = (uₖ - λₖ₋₁wₖ₋₁ - ϵₖ₋₂wₖ₋₂) / δₖ
       if iter ≥ 3
@@ -293,7 +294,7 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
         wₖ = wₖ₋₂
         @kaxpy!(n, -λₖ₋₁, wₖ₋₁, wₖ)
         @kaxpy!(n, one(FC), uₖ, wₖ)
-        @. wₖ = wₖ / δₖ
+        wₖ .= wₖ ./ δₖ
       end
 
       # Compute solution xₖ.
@@ -309,14 +310,14 @@ kwargs_usymqr = (:atol, :rtol, :itmax, :timemax, :verbose, :history, :callback,
       history && push!(AᴴrNorms, AᴴrNorm)
 
       # Compute uₖ₊₁ and uₖ₊₁.
-      @. vₖ₋₁ = vₖ # vₖ₋₁ ← vₖ
-      @. uₖ₋₁ = uₖ # uₖ₋₁ ← uₖ
+      @kcopy!(m, vₖ, vₖ₋₁)  # vₖ₋₁ ← vₖ
+      @kcopy!(n, uₖ, uₖ₋₁)  # uₖ₋₁ ← uₖ
 
       if βₖ₊₁ ≠ zero(T)
-        @. vₖ = q / βₖ₊₁ # βₖ₊₁vₖ₊₁ = q
+        vₖ .= q ./ βₖ₊₁ # βₖ₊₁vₖ₊₁ = q
       end
       if γₖ₊₁ ≠ zero(T)
-        @. uₖ = p / γₖ₊₁ # γₖ₊₁uₖ₊₁ = p
+        uₖ .= p ./ γₖ₊₁ # γₖ₊₁uₖ₊₁ = p
       end
 
       # Update directions for x.
diff --git a/test/cpu/linear_operators.jl b/test/cpu/linear_operators.jl
new file mode 100644
index 000000000..3a58c5d20
--- /dev/null
+++ b/test/cpu/linear_operators.jl
@@ -0,0 +1,16 @@
+using LinearAlgebra, SparseArrays, Test
+using Krylov, LinearOperators
+
+@testset "LinearOperators" begin
+  n = 50
+  p = 5
+
+  for T in (Float64, ComplexF64)
+    A = rand(T, n, n)
+    B = rand(T, n, p)
+    
+    opA = LinearOperator(A)
+    x, stats = block_gmres(opA, B)
+    @test stats.solved
+  end
+end
diff --git a/test/test_bicgstab.jl b/test/test_bicgstab.jl
index 6817acf3d..bde109fdd 100644
--- a/test/test_bicgstab.jl
+++ b/test/test_bicgstab.jl
@@ -62,7 +62,7 @@
       A, b, M = square_preconditioned(FC=FC)
       (x, stats) = bicgstab(A, b, M=M)
       r = b - A * x
-      resid = norm(r) / norm(b)
+      resid = norm(M * r) / norm(M * b)
       @test(resid ≤ bicgstab_tol)
       @test(stats.solved)
 
@@ -78,7 +78,7 @@
       A, b, M, N = two_preconditioners(500, 32)
       (x, stats) = bicgstab(A, b, M=M, N=N)
       r = b - A * x
-      resid = norm(r) / norm(b)
+      resid = norm(M * r) / norm(M * b)
       @test(resid ≤ bicgstab_tol)
       @test(stats.solved)
 
diff --git a/test/test_bilq.jl b/test/test_bilq.jl
index 40b9872db..672097187 100644
--- a/test/test_bilq.jl
+++ b/test/test_bilq.jl
@@ -66,11 +66,34 @@
       @test(resid ≤ bilq_tol)
       @test(stats.solved)
 
+      # Left preconditioning
+      A, b, M = square_preconditioned(FC=FC)
+      (x, stats) = bilq(A, b, M=M)
+      r = b - A * x
+      resid = norm(M * r) / norm(M * b)
+      @test(resid ≤ bilq_tol)
+      @test(stats.solved)
+
+      # Right preconditioning
+      A, b, N = square_preconditioned(FC=FC)
+      (x, stats) = bilq(A, b, N=N)
+      r = b - A * x
+      resid = norm(r) / norm(b)
+      @test(resid ≤ bilq_tol)
+      @test(stats.solved)
+
+      # Split preconditioning
+      A, b, M, N = two_preconditioners(FC=FC)
+      (x, stats) = bilq(A, b, M=M, N=N)
+      r = b - A * x
+      resid = norm(M * r) / norm(M * b)
+      @test(resid ≤ bilq_tol)
+      @test(stats.solved)
+
       # Test bᴴc == 0
       A, b, c = bc_breakdown(FC=FC)
       (x, stats) = bilq(A, b, c=c)
       @test stats.status == "Breakdown bᴴc = 0"
-
       
       # test callback function
       A, b = polar_poisson(FC=FC)
diff --git a/test/test_cgs.jl b/test/test_cgs.jl
index 832cd76c3..e4aaca4ef 100644
--- a/test/test_cgs.jl
+++ b/test/test_cgs.jl
@@ -62,7 +62,7 @@
       A, b, N = square_preconditioned(FC=FC)
       (x, stats) = cgs(A, b, N=N)
       r = b - A * x
-      resid = norm(r) / norm(b)
+      resid = norm(M * r) / norm(M * b)
       @test(resid ≤ cgs_tol)
       @test(stats.solved)
 
@@ -70,7 +70,7 @@
       A, b, M, N = two_preconditioners(FC=FC)
       (x, stats) = cgs(A, b, M=M, N=N)
       r = b - A * x
-      resid = norm(r) / norm(b)
+      resid = norm(M * r) / norm(M * b)
       @test(resid ≤ cgs_tol)
       @test(stats.solved)
 
diff --git a/test/test_qmr.jl b/test/test_qmr.jl
index 4a6b8c1c9..91fc77dd0 100644
--- a/test/test_qmr.jl
+++ b/test/test_qmr.jl
@@ -58,6 +58,30 @@
       @test(resid ≤ qmr_tol)
       @test(stats.solved)
 
+      # Left preconditioning
+      A, b, M = square_preconditioned(FC=FC)
+      (x, stats) = qmr(A, b, M=M)
+      r = b - A * x
+      resid = norm(M * r) / norm(M * b)
+      @test(resid ≤ qmr_tol)
+      @test(stats.solved)
+
+      # Right preconditioning
+      A, b, N = square_preconditioned(FC=FC)
+      (x, stats) = qmr(A, b, N=N)
+      r = b - A * x
+      resid = norm(r) / norm(b)
+      @test(resid ≤ qmr_tol)
+      @test(stats.solved)
+
+      # Split preconditioning
+      A, b, M, N = two_preconditioners(FC=FC)
+      (x, stats) = qmr(A, b, M=M, N=N)
+      r = b - A * x
+      resid = norm(M * r) / norm(M * b)
+      @test(resid ≤ qmr_tol)
+      @test(stats.solved)
+
       # Test bᴴc == 0
       A, b, c = bc_breakdown(FC=FC)
       (x, stats) = qmr(A, b, c=c)