formatting

Jutho · Dec 6, 2024 · 15b163d · 15b163d
1 parent 448acfa
commit 15b163d
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 25 deletions.
diff --git a/src/algorithms.jl b/src/algorithms.jl
@@ -314,16 +314,16 @@ struct LSMR{O<:Orthogonalizer,S<:Real} <: KrylovAlgorithm
     λ::S
     krylovdim::Int
 end
-LSMR(; orth = KrylovDefaults.orth,
-    atol = KrylovDefaults.tol,
-    btol = KrylovDefaults.tol,
-    conlim = 1/min(atol,btol),
-    maxiter = KrylovDefaults.maxiter,
-    krylovdim = KrylovDefaults.krylovdim,
-    λ = zero(atol),
-    verbosity = 0) = LSMR(orth,atol,btol,conlim,maxiter,verbosity,λ,krylovdim)
-
-
+function LSMR(; orth=KrylovDefaults.orth,
+              atol=KrylovDefaults.tol,
+              btol=KrylovDefaults.tol,
+              conlim=1 / min(atol, btol),
+              maxiter=KrylovDefaults.maxiter,
+              krylovdim=KrylovDefaults.krylovdim,
+              λ=zero(atol),
+              verbosity=0)
+    return LSMR(orth, atol, btol, conlim, maxiter, verbosity, λ, krylovdim)
+end
 
 # TODO
 """

diff --git a/src/linsolve/lsmr.jl b/src/linsolve/lsmr.jl
@@ -1,9 +1,9 @@
 # reference implementation https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/src/lsmr.jl
 function linsolve(operator, b, alg::LSMR)
-    return linsolve(operator, b, 0 * apply_adjoint(operator, b), alg)
+    return linsolve(operator, b, zerovector(apply_adjoint(operator, b)), alg)
 end;
-function linsolve(operator, b, x, alg::LSMR)
-    u = axpby!(1, b, -1, apply_normal(operator, x))
+function linsolve(operator, b, x₀, alg::LSMR)
+    u = axpby!(1, b, -1, apply_normal(operator, x₀))
     β = norm(u)
 
     # initialize GKL factorization
@@ -17,6 +17,8 @@ function linsolve(operator, b, x, alg::LSMR)
     alg.conlim > 0 ? ctol = convert(Tr, inv(alg.conlim)) : ctol = zero(Tr)
     istop = 0
 
+    x = copy(x₀)
+
     for topit in 1:(alg.maxiter)# the outermost restart loop
         # Initialize variables for 1st iteration.
         α = fact.αs[end]

diff --git a/test/linsolve.jl b/test/linsolve.jl
@@ -56,26 +56,28 @@ end
 @testset "full lsmr" begin
     @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
         @testset for orth in (cgs2, mgs2, cgsr, mgsr)
-            A = rand(T, (n,n))
-            v = rand(T,n);
-            w = rand(T,n);
-            alg = LSMR(orth = orth, krylovdim = 2*n, maxiter = 1, atol = 10*n*eps(real(T)),btol = 10*n*eps(real(T)))
-            S, info = @inferred linsolve(wrapop(A), wrapvec(v),wrapvec(w), alg)
+            A = rand(T, (n, n))
+            v = rand(T, n)
+            w = rand(T, n)
+            alg = LSMR(; orth=orth, krylovdim=2 * n, maxiter=1, atol=10 * n * eps(real(T)),
+                       btol=10 * n * eps(real(T)))
+            S, info = @inferred linsolve(wrapop(A), wrapvec(v), wrapvec(w), alg)
             @test info.converged > 0
-            @test v≈A*unwrapvec(S)+unwrapvec(info.residual)
+            @test v ≈ A * unwrapvec(S) + unwrapvec(info.residual)
         end
     end
 end
 @testset "iterative lsmr" begin
     @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
         @testset for orth in (cgs2, mgs2, cgsr, mgsr)
-            A = rand(T, (N,N))
-            v = rand(T,N);
-            w = rand(T,N);
-            alg = LSMR(orth = orth, krylovdim = N, maxiter = 50, atol = 10*N*eps(real(T)),btol = 10*N*eps(real(T)))
-            S, info = @inferred linsolve(wrapop(A), wrapvec(v),wrapvec(w), alg)
+            A = rand(T, (N, N))
+            v = rand(T, N)
+            w = rand(T, N)
+            alg = LSMR(; orth=orth, krylovdim=N, maxiter=50, atol=10 * N * eps(real(T)),
+                       btol=10 * N * eps(real(T)))
+            S, info = @inferred linsolve(wrapop(A), wrapvec(v), wrapvec(w), alg)
             @test info.converged > 0
-            @test v≈A*unwrapvec(S)+unwrapvec(info.residual)
+            @test v ≈ A * unwrapvec(S) + unwrapvec(info.residual)
         end
     end
 end