Little modifications

src/symmlq.jl

@@ -25,7 +25,8 @@ A preconditioner M may be provided in the form of a linear operator and is
 assumed to be symmetric and positive definite.
 """
 function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
-                M :: AbstractLinearOperator=opEye(), λ :: T=zero(T),
+                M :: AbstractLinearOperator=opEye(),
+                λ :: T=zero(T), transfer_to_cg :: Bool=true,
                 λest :: T=zero(T), atol :: T=√eps(T), rtol :: T=√eps(T),
                 etol :: T=√eps(T), window :: Int=0, itmax :: Int=0,
                 conlim :: T=1/√eps(T), verbose :: Bool=false) where T <: AbstractFloat
@@ -39,15 +40,15 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   MisI = isa(M, opEye)
 
   ϵM = eps(T)
-  x_lq = zeros(T, n)
+  x = zeros(T, n)
   ctol = conlim > 0 ? 1 / conlim : zero(T)
 
   # Initialize Lanczos process.
   # β₁ M v₁ = b.
   Mvold = copy(b)
   vold = M * Mvold
   β₁ = @kdot(m, vold, Mvold)
-  β₁ == 0 && return (x_lq, zeros(T, n), SimpleStats(true, true, [zero(T)], [zero(T)], "x = 0 is a zero-residual solution"))
+  β₁ == 0 && return (x, SimpleStats(true, true, [zero(T)], [zero(T)], "x = 0 is a zero-residual solution"))
   β₁ = sqrt(β₁)
   β = β₁
   @kscal!(m, one(T) / β, vold)
@@ -84,7 +85,7 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   ANorm = zero(T)
   Acond = zero(T)
   rNorm = β₁
-  rcgNorm = β * abs(ζbar)
+  rcgNorm = β₁ * abs(ζbar)
 
   xNorm = zero(T)
   xcgNorm = abs(ζbar)
@@ -138,8 +139,8 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
     # Update SYMMLQ point
     ηold = η
     ζ = ηold / γ
-    @kaxpy!(n, c * ζ, w̅, x_lq)
-    @kaxpy!(n, s * ζ, v, x_lq)
+    @kaxpy!(n, c * ζ, w̅, x)
+    @kaxpy!(n, s * ζ, v, x)
     # Update w̅
     @kaxpby!(n, -c, v, s, w̅)
 
@@ -262,8 +263,10 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   end
 
   # Compute CG point
-  @kaxpby!(m, one(T), x_lq, ζbar, w̅)
-  x_cg = w̅
+    # (xᶜ)ₖ ← (xᴸ)ₖ₋₁ + ζbarₖ * w̅ₖ
+  if solved_cg
+    @kaxpy!(m, ζbar, w̅, x)
+  end
 
   tired         && (status = "maximum number of iterations exceeded")
   ill_cond_mach && (status = "condition number seems too large for this machine")
@@ -272,5 +275,5 @@ function symmlq(A :: AbstractLinearOperator, b :: AbstractVector{T};
   zero_resid    && (status = "found approximate zero-residual solution")
 
   stats = SymmlqStats(solved, rNorms, rcgNorms, errors, errorscg, ANorm, Acond, status)
-  return (x_lq, x_cg, stats)
+  return (x, stats)
 end