Record resnorm in powm!

README.md

@@ -23,6 +23,4 @@ Pkg.add("IterativeSolvers")
 
 - [Issue #1](https://github.com/JuliaMath/IterativeSolvers.jl/issues/1) documents the implementation roadmap for iterative solvers.
 
-- The interfaces between the various algorithms are still in flux, but aim to be consistent.
-
-- [Issue #33](https://github.com/JuliaMath/IterativeSolvers.jl/issues/33) documents the implementation roadmap for randomized algorithms.
+- The interfaces between the various algorithms are still in flux, but aim to be consistent.

REQUIRE

@@ -1,4 +1,3 @@
 julia 0.6
 
 RecipesBase
-SugarBLAS

docs/make.jl

@@ -25,7 +25,6 @@ makedocs(
 			"Power method" => "eigenproblems/power_method.md",
 		],
 		"SVDL" => "svd/svdl.md",
-		"Randomized algorithms" => "randomized.md",
 		"The iterator approach" => "iterators.md",
 		"About" => [
 			"Contributing" => "about/CONTRIBUTING.md",

docs/src/index.md

@@ -4,7 +4,7 @@ IterativeSolvers.jl is a Julia package that provides efficient iterative algorit
 
 For bug reports, feature requests and questions please submit an issue. If you're interested in contributing, please see the [Contributing](@ref) guide.
 
-For more information on future methods have a look at the package [roadmap](https://github.com/JuliaLang/IterativeSolvers.jl/issues/1) on deterministic methods, for randomized algorithms check [here](https://github.com/JuliaLang/IterativeSolvers.jl/issues/33).
+For more information on future methods have a look at the package [roadmap](https://github.com/JuliaLang/IterativeSolvers.jl/issues/1).
 
 ## What method should I use for linear systems?
 
@@ -29,9 +29,3 @@ When solving **least-squares** problems we currently offer just [LSMR](@ref LSMR
 For the Singular Value Decomposition we offer [SVDL](@ref SVDL), which is the Golub-Kahan-Lanczos procedure.
 
 For eigenvalue problems we have at this point just the [Power Method](@ref PowerMethod) and some convenience wrappers to do shift-and-invert.
-
-## Randomized algorithms
-
-[Randomized algorithms](@ref Randomized) have gotten some traction lately. Some of the methods mentioned in [^Halko2011] have been implemented as well, although their quality is generally poor. Also note that many classical methods such as the subspace iteration, BiCG and recent methods like IDR(s) are also "randomized" in some sense.
-
-[^Halko2011]: Halko, Nathan, Per-Gunnar Martinsson, and Joel A. Tropp. "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions." SIAM review 53.2 (2011): 217-288.

src/IterativeSolvers.jl

@@ -5,17 +5,14 @@ eigensystems, and singular value problems.
 """
 module IterativeSolvers
 
-using SugarBLAS
-
 include("common.jl")
 include("orthogonalize.jl")
 include("history.jl")
 
-#Specialized factorizations
-include("factorization.jl")
+# Factorizations
 include("hessenberg.jl")
 
-#Linear solvers
+# Linear solvers
 include("stationary.jl")
 include("stationary_sparse.jl")
 include("cg.jl")
@@ -25,19 +22,14 @@ include("gmres.jl")
 include("chebyshev.jl")
 include("idrs.jl")
 
-#Eigensolvers
+# Eigensolvers
 include("simple.jl")
 
-#SVD solvers
+# SVD solvers
 include("svdl.jl")
 
-#Least-squares
+# Least-squares
 include("lsqr.jl")
 include("lsmr.jl")
 
-#Randomized algorithms
-include("rlinalg.jl")
-include("rsvd.jl")
-include("rsvd_fnkz.jl")
-
 end

src/bicgstabl.jl

@@ -24,8 +24,6 @@ mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, small
     M::smallMatT
 end
 
-bicgstabl_iterator(A, b, l; kwargs...) = bicgstabl_iterator!(zerox(A, b), A, b, l; initial_zero = true, kwargs...)
-
 function bicgstabl_iterator!(x, A, b, l::Int = 2;
     Pl = Identity(),
     max_mv_products = size(A, 2),
@@ -49,8 +47,7 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2;
         copy!(residual, b)
     else
         A_mul_B!(residual, A, x)
-        @blas! residual -= one(T) * b
-        @blas! residual *= -one(T)
+        residual .= b .- residual
         mv_products += 1
     end
 
@@ -91,10 +88,7 @@ function next(it::BiCGStabIterable, iteration::Int)
         β = ρ / it.σ
 
         # us[:, 1 : j] .= rs[:, 1 : j] - β * us[:, 1 : j]
-        for i = 1 : j
-            @blas! view(it.us, :, i) *= -β
-            @blas! view(it.us, :, i) += one(T) * view(it.rs, :, i)
-        end
+        it.us[:, 1 : j] .= it.rs[:, 1 : j] .- β .* it.us[:, 1 : j]
 
         # us[:, j + 1] = Pl \ (A * us[:, j])
         next_u = view(it.us, :, j + 1)
@@ -104,18 +98,15 @@ function next(it::BiCGStabIterable, iteration::Int)
         it.σ = dot(it.r_shadow, next_u)
         α = ρ / it.σ
 
-        # rs[:, 1 : j] .= rs[:, 1 : j] - α * us[:, 2 : j + 1]
-        for i = 1 : j
-            @blas! view(it.rs, :, i) -= α * view(it.us, :, i + 1)
-        end
+        it.rs[:, 1 : j] .-= α .* it.us[:, 2 : j + 1]
 
         # rs[:, j + 1] = Pl \ (A * rs[:, j])
         next_r = view(it.rs, :, j + 1)
         A_mul_B!(next_r, it.A , view(it.rs, :, j))
         A_ldiv_B!(it.Pl, next_r)
 
         # x = x + α * us[:, 1]
-        @blas! it.x += α * view(it.us, :, 1)
+        it.x .+= α .* view(it.us, :, 1)
     end
 
     # Bookkeeping

src/cg.jl

@@ -1,6 +1,6 @@
 import Base: start, next, done
 
-export cg, cg!, CGIterable, PCGIterable, cg_iterator, cg_iterator!
+export cg, cg!, CGIterable, PCGIterable, cg_iterator!, CGStateVariables
 
 mutable struct CGIterable{matT, solT, vecT, numT <: Real}
     A::matT
@@ -43,16 +43,15 @@ end
 function next(it::CGIterable, iteration::Int)
     # u := r + βu (almost an axpy)
     β = it.residual^2 / it.prev_residual^2
-    @blas! it.u *= β
-    @blas! it.u += one(eltype(it.u)) * it.r
+    it.u .= it.r .+ β .* it.u
 
     # c = A * u
     A_mul_B!(it.c, it.A, it.u)
     α = it.residual^2 / dot(it.u, it.c)
 
     # Improve solution and residual
-    @blas! it.x += α * it.u
-    @blas! it.r -= α * it.c
+    it.x .+= α .* it.u
+    it.r .-= α .* it.c
 
     it.prev_residual = it.residual
     it.residual = norm(it.r)
@@ -73,16 +72,15 @@ function next(it::PCGIterable, iteration::Int)
 
     # u := c + βu (almost an axpy)
     β = it.ρ / ρ_prev
-    @blas! it.u *= β
-    @blas! it.u += one(eltype(it.u)) * it.c
+    it.u .= it.c .+ β .* it.u
 
     # c = A * u
     A_mul_B!(it.c, it.A, it.u)
     α = it.ρ / dot(it.u, it.c)
 
     # Improve solution and residual
-    @blas! it.x += α * it.u
-    @blas! it.r -= α * it.c
+    it.x .+= α .* it.u
+    it.r .-= α .* it.c
 
     it.residual = norm(it.r)
 
@@ -92,15 +90,32 @@ end
 
 # Utility functions
 
-@inline cg_iterator(A, b, Pl = Identity(); kwargs...) = cg_iterator!(zerox(A, b), A, b, Pl; initially_zero = true, kwargs...)
+"""
+Intermediate CG state variables to be used inside cg and cg!. `u`, `r` and `c` should be of the same type as the solution of `cg` or `cg!`.
+```
+struct CGStateVariables{T,Tx<:AbstractArray{T}}
+    u::Tx
+    r::Tx
+    c::Tx
+end
+```
+"""
+struct CGStateVariables{T,Tx<:AbstractArray{T}}
+    u::Tx
+    r::Tx
+    c::Tx
+end
 
 function cg_iterator!(x, A, b, Pl = Identity();
     tol = sqrt(eps(real(eltype(b)))),
     maxiter::Int = size(A, 2),
+    statevars::CGStateVariables = CGStateVariables{eltype(x),typeof(x)}(zeros(x), similar(x), similar(x)),
     initially_zero::Bool = false
 )
-    u = zeros(x)
-    r = similar(x)
+    u = statevars.u
+    r = statevars.r
+    c = statevars.c
+    u .= zero(eltype(x))
     copy!(r, b)
 
     # Compute r with an MV-product or not.
@@ -111,8 +126,8 @@ function cg_iterator!(x, A, b, Pl = Identity();
         reltol = residual * tol # Save one dot product
     else
         mv_products = 1
-        c = A * x
-        @blas! r -= one(eltype(x)) * c
+        A_mul_B!(c, A, x)
+        r .-= c
         residual = norm(r)
         reltol = norm(b) * tol
     end
@@ -149,15 +164,16 @@ cg(A, b; kwargs...) = cg!(zerox(A, b), A, b; initially_zero = true, kwargs...)
 
 ## Keywords
 
+- `statevars::CGStateVariables`: Has 3 arrays similar to `x` to hold intermediate results;
 - `initially_zero::Bool`: If `true` assumes that `iszero(x)` so that one 
   matrix-vector product can be saved when computing the initial 
   residual vector;
 - `Pl = Identity()`: left preconditioner of the method. Should be symmetric, 
-  positive-definite like `A`.
+  positive-definite like `A`;
 - `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`;
 - `maxiter::Int = size(A,2)`: maximum number of iterations;
 - `verbose::Bool = false`: print method information;
-- `log::Bool = false`: keep track of the residual norm in each iteration;
+- `log::Bool = false`: keep track of the residual norm in each iteration.
 
 # Output
 
@@ -179,6 +195,7 @@ function cg!(x, A, b;
     tol = sqrt(eps(real(eltype(b)))),
     maxiter::Int = size(A, 2),
     log::Bool = false,
+    statevars::CGStateVariables = CGStateVariables{eltype(x), typeof(x)}(zeros(x), similar(x), similar(x)),
     verbose::Bool = false,
     Pl = Identity(),
     kwargs...
@@ -188,7 +205,7 @@ function cg!(x, A, b;
     log && reserve!(history, :resnorm, maxiter + 1)
 
     # Actually perform CG
-    iterable = cg_iterator!(x, A, b, Pl; tol = tol, maxiter = maxiter, kwargs...)
+    iterable = cg_iterator!(x, A, b, Pl; tol = tol, maxiter = maxiter, statevars = statevars, kwargs...)
     if log
         history.mvps = iterable.mv_products
     end

src/chebyshev.jl

@@ -37,26 +37,20 @@ function next(cheb::ChebyshevIterable, iteration::Int)
     else
         β = (cheb.λ_diff * cheb.α / 2) ^ 2
         cheb.α = inv(cheb.λ_avg - β)
-
-        # Almost an axpy u = c + βu
-        scale!(cheb.u, β)
-        @blas! cheb.u += one(T) * cheb.c
+        cheb.u .= cheb.c .+ β .* cheb.c
     end
 
     A_mul_B!(cheb.c, cheb.A, cheb.u)
     cheb.mv_products += 1
 
-    @blas! cheb.x += cheb.α * cheb.u
-    @blas! cheb.r -= cheb.α * cheb.c
+    cheb.x .+= cheb.α .* cheb.u
+    cheb.r .-= cheb.α .* cheb.c
 
     cheb.resnorm = norm(cheb.r)
 
     cheb.resnorm, iteration + 1
 end
 
-chebyshev_iterable(A, b, λmin::Real, λmax::Real; kwargs...) =
-    chebyshev_iterable!(zerox(A, b), A, b, λmin, λmax; kwargs...)
-
 function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
     tol = sqrt(eps(real(eltype(b)))), maxiter = size(A, 2), Pl = Identity(), initially_zero = false)
 
@@ -76,7 +70,7 @@ function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
         mv_products = 0
     else
         A_mul_B!(c, A, x)
-        @blas! r -= one(T) * c
+        r .-= c
         resnorm = norm(r)
         reltol = tol * norm(b)
         mv_products = 1