diff --git a/README.md b/README.md index feec0f06..cdde6710 100644 --- a/README.md +++ b/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. \ No newline at end of file diff --git a/REQUIRE b/REQUIRE index ca488e5f..0284cb21 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,3 @@ julia 0.6 RecipesBase -SugarBLAS diff --git a/docs/make.jl b/docs/make.jl index 6254d794..422b3487 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/index.md b/docs/src/index.md index 0b92897e..2881c612 100644 --- a/docs/src/index.md +++ b/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. \ No newline at end of file diff --git a/docs/src/randomized.md b/docs/src/randomized.md deleted file mode 100644 index 707b84a1..00000000 --- a/docs/src/randomized.md +++ /dev/null @@ -1,35 +0,0 @@ -# [Randomized algorithms](@id Randomized) - -The methods below are based on [^Halko2011]. - -```@docs -reig -rsvdfact -rsvd_fnkz -``` - -## Condition number estimate - -```@docs -rcond -``` - -## Extremal eigenvalue estimates - -```@docs -reigmin -reigmax -``` - -## Norm estimate - -```@docs -rnorm -rnorms -``` - -[^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. - -[^Dixon1983]: Dixon, John D. "Estimating extremal eigenvalues and condition numbers of matrices." SIAM Journal on Numerical Analysis 20.4 (1983): 812-814. - -[^Liberty2007]: Liberty, Edo, et al. "Randomized algorithms for the low-rank approximation of matrices." Proceedings of the National Academy of Sciences 104.51 (2007): 20167-20172. \ No newline at end of file diff --git a/src/IterativeSolvers.jl b/src/IterativeSolvers.jl index 78c3b51e..07696980 100644 --- a/src/IterativeSolvers.jl +++ b/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 diff --git a/src/bicgstabl.jl b/src/bicgstabl.jl index 8db48499..7f1a2302 100644 --- a/src/bicgstabl.jl +++ b/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,10 +98,7 @@ 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) @@ -115,7 +106,7 @@ function next(it::BiCGStabIterable, iteration::Int) 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 diff --git a/src/cg.jl b/src/cg.jl index 14c91d40..26db481f 100644 --- a/src/cg.jl +++ b/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 diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 3fd5d58d..dcb328cc 100644 --- a/src/chebyshev.jl +++ b/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 diff --git a/src/factorization.jl b/src/factorization.jl deleted file mode 100644 index c370915e..00000000 --- a/src/factorization.jl +++ /dev/null @@ -1,90 +0,0 @@ -############################ -# Specialized factorizations -############################ - -export idfact - -""" -An interpolative decomposition. - -For a matrix `A`, the interpolative decomposition `F` contains the matrices `B` -and `P` computed by `idfact()`. See the documentation of `idfact()` for details. - -# References - -\\cite{Cheng2005, Liberty2007} -""" -struct Interpolative{T} <: Factorization{T} - B :: AbstractMatrix{T} - P :: AbstractMatrix{T} -end - -""" - idfact(A, k, l) - -Compute and return the interpolative decomposition of `A`: A ≈ B * P - -Where: -* `B`'s columns are a subset of the columns of `A` -* some subset of `P`'s columns are the `k x k` identity, no entry of `P` exceeds magnitude 2, and -* ||B * P - A|| ≲ σ(A, k+1), the (`k+1`)st singular value of `A`. - -# Arguments - -`A`: Matrix to factorize - -`k::Int`: Number of columns of A to return in B - -`l::Int`: Length of random vectors to project onto - -# Output - -`(::Interpolative)`: interpolative decomposition. - -# Implementation note - -This is a hacky version of the algorithms described in \\cite{Liberty2007} -and \\cite{Cheng2005}. The former refers to the factorization (3.1) of the -latter. However, it is not actually necessary to compute this -factorization in its entirely to compute an interpolative decomposition. - -Instead, it suffices to find some permutation of the first k columns of Y = -R * A, extract the subset of A into B, then compute the P matrix as B\\A -which will automatically compute P using a suitable least-squares -algorithm. - -The approximation we use here is to compute the column pivots of Y, -rather then use the true column pivots as would be computed by a column- -pivoted QR process. - -# References - -\\cite[Algorithm I]{Liberty2007} - -```bibtex -@article{Cheng2005, - author = {Cheng, H and Gimbutas, Z and Martinsson, P G and Rokhlin, V}, - doi = {10.1137/030602678}, - issn = {1064-8275}, - journal = {SIAM Journal on Scientific Computing}, - month = jan, - number = {4}, - pages = {1389--1404}, - title = {On the Compression of Low Rank Matrices}, - volume = {26}, - year = {2005} -} -``` -""" -function idfact(A, k::Int, l::Int) - m, n = size(A) - R = randn(l, m) - Y = R * A #size l x n - - #Compute column pivots of first k columns of Y - maxvals = map(j->maximum(abs(view(Y, :, j))), 1:n) - piv = sortperm(maxvals, rev=true)[1:k] - - B = A[:, piv] - Interpolative(B, B\A) -end diff --git a/src/gmres.jl b/src/gmres.jl index ce231bc6..7a19b8fc 100644 --- a/src/gmres.jl +++ b/src/gmres.jl @@ -2,7 +2,7 @@ import Base: start, next, done export gmres, gmres! -mutable struct ArnoldiDecomp{T, matT} +struct ArnoldiDecomp{T, matT} A::matT V::Matrix{T} # Orthonormal basis vectors H::Matrix{T} # Hessenberg matrix @@ -98,8 +98,6 @@ function next(g::GMRESIterable, iteration::Int) g.residual.current, iteration + 1 end -gmres_iterable(A, b; kwargs...) = gmres_iterable!(zerox(A, b), A, b; initially_zero = true, kwargs...) - function gmres_iterable!(x, A, b; Pl = Identity(), Pr = Identity(), @@ -223,14 +221,14 @@ function init!(arnoldi::ArnoldiDecomp{T}, x, b, Pl, Ax; initially_zero::Bool = f # Potentially save one MV product if !initially_zero A_mul_B!(Ax, arnoldi.A, x) - @blas! first_col -= one(T) * Ax + first_col .-= Ax end A_ldiv_B!(Pl, first_col) # Normalize β = norm(first_col) - @blas! first_col *= one(T) / β + first_col .*= inv(β) β end @@ -244,7 +242,7 @@ function solve_least_squares!(arnoldi::ArnoldiDecomp{T}, β, k::Int) where {T} rhs = zeros(T, k) rhs[1] = β - H = Hessenberg(view(arnoldi.H, 1 : k, 1 : k - 1)) + H = FastHessenberg(view(arnoldi.H, 1 : k, 1 : k - 1)) A_ldiv_B!(H, rhs) rhs @@ -261,7 +259,7 @@ function update_solution!(x, y, arnoldi::ArnoldiDecomp{T}, Pr, k::Int, Ax) where # Computing x ← x + Pr \ (V * y) and use Ax as a work space A_mul_B!(Ax, view(arnoldi.V, :, 1 : k - 1), y) A_ldiv_B!(Pr, Ax) - @blas! x += one(T) * Ax + x .+= Ax end function expand!(arnoldi::ArnoldiDecomp, Pl::Identity, Pr::Identity, k::Int, Ax) diff --git a/src/hessenberg.jl b/src/hessenberg.jl index 0ce06eea..f1f11dc3 100644 --- a/src/hessenberg.jl +++ b/src/hessenberg.jl @@ -1,18 +1,18 @@ import Base.LinAlg: Givens, givensAlgorithm import Base.A_ldiv_B! -mutable struct Hessenberg{T<:AbstractMatrix} +struct FastHessenberg{T<:AbstractMatrix} H::T # H is assumed to be Hessenberg of size (m + 1) × m end -@inline Base.size(H::Hessenberg, args...) = size(H.H, args...) +@inline Base.size(H::FastHessenberg, args...) = size(H.H, args...) """ Solve Hy = rhs for a non-square Hessenberg matrix. Note that `H` is also modified as is it converted to an upper triangular matrix via Given's rotations """ -function A_ldiv_B!(H::Hessenberg, rhs) +function A_ldiv_B!(H::FastHessenberg, rhs) # Implicitly computes H = QR via Given's rotations # and then computes the least-squares solution y to # |Hy - rhs| = |QRy - rhs| = |Ry - Q'rhs| diff --git a/src/history.jl b/src/history.jl index 70b79a5c..758e5922 100644 --- a/src/history.jl +++ b/src/history.jl @@ -1,6 +1,6 @@ using RecipesBase -import Base: getindex, setindex!, push!, keys +import Base: getindex, setindex!, push!, keys, show export ConvergenceHistory export nprods, niters, nrests @@ -90,6 +90,11 @@ const RestartedHistory{T} = ConvergenceHistory{T, Int} # Functions # ############# +function show(io::IO, ch::ConvergenceHistory) + print(io, ch.isconverged ? "Converged" : "Not converged", + " after ", ch.iters, " iterations.") +end + """ getindex(ch, s) diff --git a/src/idrs.jl b/src/idrs.jl index 9338c5a4..21d78fc1 100644 --- a/src/idrs.jl +++ b/src/idrs.jl @@ -45,7 +45,7 @@ function idrs!(x, A, b; history = ConvergenceHistory(partial=!log) history[:tol] = tol reserve!(history,:resnorm, maxiter) - idrs_method!(history, x, linsys_op, (A,), b, s, tol, maxiter; kwargs...) + idrs_method!(history, x, A, b, s, tol, maxiter; kwargs...) log && shrink!(history) log ? (x, history) : x end @@ -67,14 +67,12 @@ end om end -@inline linsys_op(x, A) = A*x - -function idrs_method!(log::ConvergenceHistory, X, op, args, C::T, +function idrs_method!(log::ConvergenceHistory, X, A, C::T, s::Number, tol::Number, maxiter::Number; smoothing::Bool=false, verbose::Bool=false ) where {T} verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm") - R = C - op(X, args...)::T + R = C - A*X normR = vecnorm(R) iter = 1 @@ -111,32 +109,26 @@ function idrs_method!(log::ConvergenceHistory, X, op, args, C::T, # Solve small system and make v orthogonal to P c = LowerTriangular(M[k:s,k:s])\f[k:s] - @blas! V = G[k] - @blas! V *= c[1] + V .= c[1] .* G[k] + Q .= c[1] .* U[k] - @blas! Q = U[k] - @blas! Q *= c[1] for i = k+1:s - @blas! V += c[i-k+1]*G[i] - @blas! Q += c[i-k+1]*U[i] + V .+= c[i-k+1] .* G[i] + Q .+= c[i-k+1] .* U[i] end # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j + V .= R .- V - #V = R - V - @blas! V *= -1. - @blas! V += R - - @blas! U[k] = Q - @blas! U[k] += om*V - G[k] = op(U[k], args...) + U[k] .= Q .+ om .* V + A_mul_B!(G[k], A, U[k]) # Bi-orthogonalise the new basis vectors for i in 1:k-1 alpha = vecdot(P[i],G[k])/M[i,i] - @blas! G[k] += -alpha*G[i] - @blas! U[k] += -alpha*U[i] + G[k] .-= alpha .* G[i] + U[k] .-= alpha .* U[i] end # New column of M = P'*G (first k-1 entries are zero) @@ -148,22 +140,17 @@ function idrs_method!(log::ConvergenceHistory, X, op, args, C::T, # Make r orthogonal to q_i, i = 1..k beta = f[k]/M[k,k] - @blas! R += -beta*G[k] - @blas! X += beta*U[k] + R .-= beta .* G[k] + X .+= beta .* U[k] normR = vecnorm(R) if smoothing - # T_s = R_s - R - @blas! T_s = R_s - @blas! T_s += (-1.)*R + T_s .= R_s .- R gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s) - @blas! R_s += -gamma*T_s - # X_s = X_s - gamma*(X_s - X) - @blas! T_s = X_s - @blas! T_s += (-1.)*X - @blas! X_s += -gamma*T_s + R_s .-= gamma .* T_s + X_s .-= gamma .* (X_s .- X) normR = vecnorm(R_s) end @@ -175,41 +162,36 @@ function idrs_method!(log::ConvergenceHistory, X, op, args, C::T, return X end if k < s - f[k+1:s] = f[k+1:s] - beta*M[k+1:s,k] + f[k+1:s] .-= beta*M[k+1:s,k] end iter += 1 end # Now we have sufficient vectors in G_j to compute residual in G_j+1 # Note: r is already perpendicular to P so v = r - @blas! V = R - Q = op(V, args...)::T + copy!(V, R) + A_mul_B!(Q, A, V) om = omega(Q, R) - @blas! R += -om*Q - @blas! X += om*V + R .-= om .* Q + X .+= om .* V normR = vecnorm(R) if smoothing - # T_s = R_s - R - @blas! T_s = R_s - @blas! T_s += (-1.)*R + T_s .= R_s .- R gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s) - @blas! R_s += -gamma*T_s - # X_s = X_s - gamma*(X_s - X) - @blas! T_s = X_s - @blas! T_s += (-1.)*X - @blas! X_s += -gamma*T_s + R_s .-= gamma .* T_s + X_s .-= gamma .* (X_s .- X) normR = vecnorm(R_s) end iter += 1 - nextiter!(log,mvps=1) + nextiter!(log, mvps=1) push!(log, :resnorm, normR) end if smoothing - @blas! X = X_s + copy!(X, X_s) end verbose && @printf("\n") setconv(log, 0<=normR 0 && @blas! u *= inv(β) + u .*= inv(β) Ac_mul_B!(1, A, u, 0, v) α = norm(v) - α > 0 && @blas! v *= inv(α) + v .*= inv(α) log[:atol] = atol log[:btol] = btol @@ -126,7 +126,7 @@ function lsmr_method!(log::ConvergenceHistory, x, A, b, v, h, hbar; cbar = one(Tr) sbar = zero(Tr) - @blas! h = v + copy!(h, v) fill!(hbar, zero(Tr)) # Initialize variables for estimation of ||r||. @@ -162,10 +162,10 @@ function lsmr_method!(log::ConvergenceHistory, x, A, b, v, h, hbar; β = norm(u) if β > 0 log.mtvps+=1 - @blas! u *= inv(β) + u .*= inv(β) Ac_mul_B!(1, A, u, -β, v) α = norm(v) - α > 0 && @blas! v *= inv(α) + v .*= inv(α) end # Construct rotation Qhat_{k,2k+1}. @@ -193,11 +193,9 @@ function lsmr_method!(log::ConvergenceHistory, x, A, b, v, h, hbar; ζbar = - sbar * ζbar # Update h, h_hat, x. - @blas! hbar *= - θbar * ρ / (ρold * ρbarold) - @blas! hbar += h - @blas! x += (ζ / (ρ * ρbar))*hbar - @blas! h *= - θnew / ρ - @blas! h += v + hbar .= hbar .* (-θbar * ρ / (ρold * ρbarold)) .+ h + x .+= (ζ / (ρ * ρbar)) * hbar + h .= h .* (-θnew / ρ) .+ v ############################################################################## ## diff --git a/src/lsqr.jl b/src/lsqr.jl index 0a9f1e33..e2af672f 100644 --- a/src/lsqr.jl +++ b/src/lsqr.jl @@ -125,12 +125,12 @@ function lsqr_method!(log::ConvergenceHistory, x, A, b; alpha = zero(Tr) if beta > 0 log.mtvps=1 - @blas! u *= inv(beta) + u .*= inv(beta) Ac_mul_B!(v,A,u) alpha = norm(v) end if alpha > 0 - @blas! v *= inv(alpha) + v .*= inv(alpha) end w = copy(v) wrho = similar(w) @@ -158,21 +158,19 @@ function lsqr_method!(log::ConvergenceHistory, x, A, b; # Note that the following three lines are a band aid for a GEMM: X: C := αAB + βC. # This is already supported in A_mul_B! for sparse and distributed matrices, but not yet dense A_mul_B!(tmpm, A, v) - @blas! u *= -alpha - @blas! u += one(eltype(tmpm))*tmpm + u .= -alpha .* u .+ tmpm beta = norm(u) if beta > 0 log.mtvps+=1 - @blas! u *= inv(beta) + u .*= inv(beta) Anorm = sqrt(abs2(Anorm) + abs2(alpha) + abs2(beta) + dampsq) # Note that the following three lines are a band aid for a GEMM: X: C := αA'B + βC. # This is already supported in Ac_mul_B! for sparse and distributed matrices, but not yet dense Ac_mul_B!(tmpn, A, u) - @blas! v *= -beta - @blas! v += one(eltype(tmpn))*tmpn + v .= -beta .* v .+ tmpn alpha = norm(v) if alpha > 0 - @blas! v *= inv(alpha) + v .*= inv(alpha) end end @@ -199,11 +197,9 @@ function lsqr_method!(log::ConvergenceHistory, x, A, b; t1 = phi /rho t2 = - theta/rho - @blas! x += t1*w - @blas! w *= t2 - @blas! w += one(t2)*v - @blas! wrho = w - @blas! wrho *= inv(rho) + x .+= t1*w + w = t2 .* w .+ v + wrho .= w .* inv(rho) ddnorm += norm(wrho) # Use a plane rotation on the right to eliminate the diff --git a/src/minres.jl b/src/minres.jl index e6df12f0..8a72a65a 100644 --- a/src/minres.jl +++ b/src/minres.jl @@ -36,8 +36,6 @@ mutable struct MINRESIterable{matT, solT, vecT <: DenseVector, smallVecT <: Dens resnorm::realT end -minres_iterable(A, b; kwargs...) = minres_iterable!(zerox(A, b), A, b; initially_zero = true, kwargs...) - function minres_iterable!(x, A, b; initially_zero::Bool = false, skew_hermitian::Bool = false, diff --git a/src/orthogonalize.jl b/src/orthogonalize.jl index 001db757..73e9bb6e 100644 --- a/src/orthogonalize.jl +++ b/src/orthogonalize.jl @@ -27,12 +27,12 @@ function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, projection_size = norm(correction) # w = w - V * correction BLAS.gemv!('N', -one(T), V, correction, one(T), w) - @blas! h += one(T) * correction + h .+= correction nrm = norm(w) end # Normalize; note that we already have norm(w). - scale!(w, one(T) / nrm) + w .*= inv(nrm) nrm end @@ -44,7 +44,7 @@ function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, nrm = norm(w) # Normalize - scale!(w, one(T) / nrm) + w .*= inv(nrm) nrm end @@ -53,12 +53,12 @@ function orthogonalize_and_normalize!(V::StridedVector{Vector{T}}, w::StridedVec # Orthogonalize using BLAS-1 ops for i = 1 : length(V) h[i] = dot(V[i], w) - @blas! w -= h[i] * V[i] + w .-= h[i] .* V[i] end # Normalize nrm = norm(w) - scale!(w, one(T) / nrm) + w .*= inv(nrm) nrm end @@ -68,11 +68,11 @@ function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, for i = 1 : size(V, 2) column = view(V, :, i) h[i] = dot(column, w) - @blas! w -= h[i] * column + w .-= h[i] .* column end nrm = norm(w) - scale!(w, one(T) / nrm) + w .*= inv(nrm) nrm end \ No newline at end of file diff --git a/src/rlinalg.jl b/src/rlinalg.jl deleted file mode 100644 index fe64eeb5..00000000 --- a/src/rlinalg.jl +++ /dev/null @@ -1,314 +0,0 @@ -################################################################# -# Randomized estimators of elementary linear algebraic quantities -################################################################# - -export rcond, reigmax, reigmin, rnorm, rnorms - -import Base: * - - -""" - randnn(el, m) - randnn(el, m, n) - -Compute randomized gaussian matrix normalized by column. - -# Arguments - -- `el::Type`: element type. -- `m::Int`: number of rows. -- `n::Int`: number of columns. - -## Keywords - -- `normalize::Bool = true`: normalize output. - -# Output - -**Without `n`** - -- `Ω`: vector containing Gaussian random numbers of type `el`. - -**With `n`** - -- `Ω`: matrix of dimensions `m` x `n` containing Gaussian random numbers of - type `el`. -""" -function randnn(el, m::Int, normalize::Bool=true) - if el <: Real - Ω = randn(m) - elseif el <: Complex - Ω = randn(m) + im*randn(m) - else - throw(ArgumentError("Unsupported element type: $el")) - end - normalize ? Ω/norm(Ω) : Ω -end -function randnn(el, m::Int, n::Int, normalize::Bool=true) - if el <: Real - Ω = randn(m, n) - elseif el <: Complex - Ω = randn(m, n) + im*randn(m, n) - else - throw(ArgumentError("Unsupported element type: $el")) - end - normalize || return Ω - for i=1:n - Ω[:, i] /= norm(view(Ω, :, i)) - end - Ω -end - -""" - rnorm(A, mvps) - -Compute a probabilistic upper bound on the norm of a matrix `A`. -`‖A‖ ≤ α √(2/π) maxᵢ ‖Aωᵢ‖` with probability `p=α^(-mvps)`. - -# Arguments - -- `A`: matrix whose norm to estimate. -- `mvps::Int`: number of matrix-vector products to compute. - -## Keywords - -- `p::Real=0.05`: probability of upper bound failing. - -# Output - -Estimate of ‖A‖. - -See also [`rnorms`](@ref) for a different estimator that uses -premultiplying by both `A` and `A'`. - -# References - -Lemma 4.1 of Halko2011 -""" -function rnorm(A, r::Int, p::Real=0.05) - @assert 0 m - throw(ArgumentError("Cannot find $l linearly independent vectors of $m x $n matrix")) - end - Ω = randn(n, l) - Y = A*Ω - Q = full(qrfact!(Y)[:Q]) - @assert m==size(Q, 1) - @assert l==size(Q, 2) - return Q -end - -""" - rrange_adaptive(A, r, ϵ=eps(); maxiter=10) - -Compute an orthogonal basis for a subspace of `A` of dimension `l` using adaptive -randomized rangefinding. - -Similar to `rrange`, but determines the oversampling parameter adaptively given -a threshold `ϵ`. - -# Arguments - -- `A`: input matrix. Must support `size(A)` and premultiply. -- `r::Integer`: number of basis vectors to compute. -- `ϵ::Real = eps()`: threshold to determine adaptive fitting. - -## Keywords - -- `maxiter::Int = 10`: maximum number of iterations to run. - -# Output - -`::Matrix`: matrix of dimension `size(A,1)` x l containing the basis -vectors of the computed subspace of `A`. - -# References - -Algorithm 4.2 of [^Halko2011] -""" -function rrange_adaptive(A, r::Integer, ϵ::Real=eps(); maxiter::Int=10) - m, n = size(A) - - Ω = randn(n,r) - Y = A*Ω - Q = zeros(m,0) - - const tol=ϵ/(10*√(2/π)) - for j=1:maxiter - Tol = maximum([norm(Y[:,i]) for i=j:(j+r-1)]) - Tol > tol || break - - y = view(Y,:,j) - y = Y[:,j] = y - Q*(Q'y) - q = y/norm(y) - Q = [Q q] - - ω = randn(n) - y = A*ω - y = y - Q*(Q'y) - Y = [Y y] - Yb = view(Y, :, (j+1):(j+r-1)) - Yb = Y[:, (j+1):(j+r-1)] = Yb - q * q'Yb - - j==maxiter && warn("Maximum number of iterations reached with norm $Tol > $tol") - end - Q -end - -""" - rrange_si(A, l; At=A', q=0) - -Compute an orthonormal basis for a subspace of `A` of dimension `l` using -randomized rangefinding by subspace iteration. - -*Note:* Running with `q=0` is functionally equivalent to `rrange`. - -# Arguments - -- `A`: input matrix. Must support `size(A)` and premultiply. -- `l::Int`: number of basis vectors to compute. - -## Keywords - -- `At = A'`: transpose of `A`. -- `q::Int = 0`: number of subspace iterations. -- `p` : oversampling parameter. The number of extra basis vectors to use - in the computation, which get discarded at the end. - -# Output - -- `::Matrix`: A dense matrix of dimension `size(A,1) x l` containing the basis - vectors of the computed subspace of `A`. - -# Implementation note - -Whereas the Reference recommends classical Gram-Schmidt with double -reorthogonalization, we instead compute the basis with `qrfact()`, which -for dense A computes the QR factorization using Householder reflectors. - -# References - -Algorithm 4.4 of [^Halko2011] -""" -function rrange_si(A, l::Int; At=A', q::Int=0) - basis=x->full(qrfact(x)[:Q]) - n = size(A, 2) - Ω = randn(n,l+p) - Y = A*Ω - Q = basis(Y) - for iter=1:q #Do some power iterations - Ỹ=At*Q - Q̃=basis(Ỹ) - Y=A*Q̃ - Q=basis(Y) - end - Q -end - -""" - rrange_f(A, l) - -Compute an orthonormal basis for a subspace of `A` of dimension `l` using -naive randomized rangefinding using stochastic randomized Fourier transforms. - -*Note:* similar to `rrange`, but does not use gaussian random matrices. - -# Arguments - -- `A` : input matrix. Must support `size(A)` and premultiply. -- `l::Int` : number of basis vectors to compute. -- `p::Int = 0` : oversampling parameter. The number of extra basis vectors to use - in the computation, which get discarded at the end. - -# Output - -- `::Matrix`: matrix of dimension `size(A,1)` x `l` containing the basis - vectors of the computed subspace of `A`. - -# Implementation note - -Whereas the Reference recommends classical Gram-Schmidt with double -reorthogonalization, we instead compute the basis with `qrfact()`, which -for dense `A` computes the QR factorization using Householder reflectors. - -# References - -Algorithm 4.5 of [^Halko2011] -""" -function rrange_f(A, l::Int) - n = size(A, 2) - Ω = srft(l+p) - Y = A*Ω - Q = full(qrfact!(Y)[:Q]) -end - -""" - svdfact_restricted(A, Q, n) - -Compute the SVD factorization of `A` restricted to the subspace spanned by `Q` -using exact projection. - -# Arguments - -- `A`: input matrix. Must support postmultiply. -- `Q`: matrix containing basis vectors of the subspace whose restriction to is - desired. - -# Output - -`::SVD`: singular value decomposition. - -# References - -Algorithm 5.1 of [^Halko2011] -""" -function svdfact_restricted(A, Q, n::Int) - B=Q'A - S=svdfact!(B) - SVD((Q*S[:U])[:, 1:n], S[:S][1:n], S[:Vt][1:n, :]) -end - -""" - svdvals_restricted(A, Q, n) - -Compute the singular values of `A` restricted to the subspace spanned by `Q` -using exact projection. - -# Arguments - -- `A`: input matrix. Must support postmultiply. -- `Q`: matrix containing basis vectors of the subspace whose restriction to is - desired. - -# Output - -- `::Vector`: estimated singular values of `A`. - -# References - -Algorithm 5.1 of [^Halko201] -""" -function svdvals_restricted(A, Q, n::Int) - B=Q'A - S=svdvals!(B)[1:n] -end - -""" - svdfact_re(A, Q) - -Compute the SVD factorization of `A` restricted to the subspace spanned by `Q` -using row extraction. - -*Note:* Remark 5.2 of [^Halko2011] recommends input of `Q` of the form `Q=A*Ω` -where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`. - -# Arguments - -- `A`: input matrix. Must support postmultiply -- `Q`: matrix containing basis vectors of the subspace whose restriction to is - desired. Need not be orthogonal or normalized. - -# Output - -- `::SVD`: singular value decomposition. - -# See also - -A faster but less accurate variant of [`svdfact_restricted`](@ref) which uses the -interpolative decomposition `idfact`. - -# References - -Algorithm 5.2 of [^Halko2011] -""" -function svdfact_re(A, Q) - F = idfact(Q) - X, J = F[:B], F[:P] - R′, W′= qr(A[J, :]) - Z = X*R′ - S=svdfact(Z) - SVD(S[:U], S[:S], S[:Vt]*W′) -end - -""" - eigfact_restricted(A, Q) - -Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace -spanned by `Q` using row extraction. - -# Arguments - -- `A::Hermitian`: input matrix. Must be `Hermitian` and support pre- and post-multiply. -- `Q`: orthonormal matrix containing basis vectors of the subspace whose - restriction to is desired. - -# Output - -- `::Base.LinAlg.Eigen`: eigen factorization. - -# References - -Algorithm 5.3 of [^Halko2011] -""" -function eigfact_restricted(A::Hermitian, Q) - B = Q'A*Q - E = eigfact!(B) - Eigen(E[:values], Q*E[:vectors]) -end - -""" - eigfact_re(A, Q) - -Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace -spanned by `Q` using row extraction. - -*Note:* Remark 5.2 of [^Halko2011] recommends input of `Q` of the form `Q=A*Ω` -where `Ω` is a sample computed by `randn(n,l)` or even `srft(l)`. - -# Arguments - -- `A::Hermitian`: input matrix. Must be `Hermitian` and support pre- and post-multiply. -- `Q`: matrix containing basis vectors of the subspace whose restriction to is - desired. Need not be orthogonal or normalized. - -# Output - -- `::Base.LinAlg.Eigen`: eigen factorization. - -# See also - -A faster but less accurate variant of `eigfact_restricted()` which uses the -interpolative decomposition `idfact()`. - -# References - -Algorithm 5.4 of [^Halko2011] -""" -function eigfact_re(A::Hermitian, Q) - X, J = idfact(Q) - F = qrfact!(X) - V, R = F[:Q], F[:R] - Z=R*A[J, J]*R' - E=eigfact(Z) - Eigen(E[:values], V*E[:vectors]) -end - -""" - eigfact_nystrom(A, Q) - -Compute the spectral (`Eigen`) factorization of `A` restricted to the subspace -spanned by `Q` using the Nyström method. - -# Arguments - -- `A`: input matrix. Must be positive semidefinite. -- `Q`: orthonormal matrix containing basis vectors of the subspace whose - restriction to is desired. - -# Output - -- `::Base.LinAlg.Eigen`: eigen factorization. - -# See also - -More accurate than [`eigfact_restricted`](@ref) but is restricted to matrices -that can be Cholesky decomposed. - -# References - -Algorithm 5.5 of [^Halko2011] -""" -function eigfact_nystrom(A, Q) - B₁=A*Q - B₂=Q'*B₁ - C=cholfact!(B₂) - F=B₁*inv(C) - S=svdfact!(F) - Eigen(S[:S].^2, S[:U]) -end - -""" - eigfact_onepass(A, Ω) - -Compute the spectral (`Eigen`) factorization of `A` using only one matrix -product involving `A`. - -# Arguments - -- `A::Hermitian`: input matrix. -- `Ω`: sample matrix for the column space, e.g. `randn(n, l)` or `srft(l)`. -- `Ω̃;`: sample matrix for the row space. Not neeeded for `Hermitian` matrices. -- `At = A'`: computes transpose of input matrix. - -# Output - -- `::Base.LingAlg.Eigen`: eigen factorization. - -# References - -Algorithm 5.6 of [^Halko2011] -""" -function eigfact_onepass(A::Hermitian, Ω) - Y=A*Ω; Q = full(qrfact!(Y)[:Q]) - B=(Q'Y)\(Q'Ω) - E=eigfact!(B) - Eigen(E[:values], Q*E[:vectors]) -end -function eigfact_onepass(A, Ω, Ω̃; At=A') - Y=A *Ω; Q = full(qrfact!(Y)[:Q]) - Ỹ=At*Ω; Q̃ = full(qrfact!(Ỹ)[:Q]) - #Want least-squares solution to (5.14 and 5.15) - B=(Q'Y)\(Q̃'Ω) - B̃=(Q̃'Ỹ)\(Q'Ω̃) - #Here is a very very very hacky way to solve the problem - B=0.5(B + B̃') - E=eigfact!(B) - Eigen(E[:values], Q*E[:vectors]) -end - -""" - reig(A, l) - -Compute the spectral (`Eigen`) decomposition of `A` using a randomized -algorithm. - -# Arguments - -- `A`: input matrix. -- `l::Int`: number of eigenpairs to find. - -# Output - -- `::Base.LinAlg.Eigen`: eigen decomposition. - -# Implementation note - -This is a wrapper around `eigfact_onepass()` which uses the randomized -samples found using `srft(l)`. -""" -reig(A::Hermitian, l::Int) = eigfact_onepass(A, srft(l)) -reig(A, l::Int) = eigfact_onepass(A, srft(l), srft(l)) diff --git a/src/rsvd_fnkz.jl b/src/rsvd_fnkz.jl deleted file mode 100644 index 2b6f9881..00000000 --- a/src/rsvd_fnkz.jl +++ /dev/null @@ -1,101 +0,0 @@ - -export rsvd_fnkz - -struct OuterProduct{T} - X :: Matrix{T} - Y :: Matrix{T} -end - -X ⊗ Y = OuterProduct{promote_type(eltype(X), eltype(Y))}(X, Y) - -""" - rsvd_fnkz(A, k) - -Compute the randomized SVD by iterative refinement from randomly selected -columns/rows [^Friedland2006]. - -# Arguments - -- `A`: matrix whose SVD is desired; -- `k::Int`: desired rank of approximation (`k ≤ min(m, n)`). - -## Keywords - -- `l::Int = k`: number of columns/rows to sample at each iteration (`1 ≤ l ≤ k`); -- `N::Int = minimum(size(A))`: maximum number of iterations; -- `ϵ::Real = prod(size(A))*eps()`: relative threshold for convergence, as - measured by growth of the spectral norm; -- `method::Symbol = :eig`: problem to solve. - 1. `:eig`: eigenproblem. - 2. `:svd`: singular problem. -- `verbose::Bool = false`: print convergence information at each iteration. - -# Return value - -SVD object of `rank ≤ k`. - -[^Friedland2006]: - Friedland, Shmuel, et al. "Fast Monte-Carlo low rank approximations for - matrices." System of Systems Engineering, 2006 IEEE/SMC International - Conference on. IEEE, 2006. -""" -function rsvd_fnkz(A, k::Int; - l::Int=k, N::Int=minimum(size(A)), verbose::Bool=false, - method::Symbol=:eig, - ϵ::Real=prod(size(A))*eps(real(float(one(eltype(A)))))) - - const m, n = size(A) - const dosvd = method == :svd - @assert 1 ≤ l ≤ k ≤ min(m, n) - - #Initialize k-rank approximation B₀ according to (III.9) - const tallandskinny = m ≥ n - B₀ = if tallandskinny - #Select k columns of A - π = randperm(n)[1:k] - A[:, π] - else - warn("Not implemented") - #Select k rows of A - π = randperm(m)[1:k] - A[π, :] - end - X, RRR = qr(B₀) - X = X[:, abs.(diag(RRR)) .> ϵ] #Remove linear dependent columns - B = tallandskinny ? X ⊗ A'X : X ⊗ X'A - - oldnrmB = 0.0 - Λ = 0 - X̃ = 0 - V = 0 - for t=1:N - π = randperm(k)[1:l] - #Update B using Theorem 2.4 - X, RRR = qr([B.X A[:, π]]) #We are doing more work than needed here - X = X[:, abs.(diag(RRR)) .> ϵ] #Remove linearly dependent columns - Y = A'X - if dosvd - S = svdfact!(Y) - Λ = S[:S].^2 - B = (X * S[:V]) ⊗ (S[:U]*Diagonal(S[:S])) - else - E = eigfact!(Symmetric(Y'Y)) - π = sortperm(E[:values], rev=true)[1:min(length(E[:values]), k)] - Λ = E[:values][π] - O = E[:vectors][:, π] #Eq. 2.6 - X̃ = X * O - B = X̃ ⊗ A'X̃ - end - - length(Λ)==0 && break - - nrmB = √Λ[1] - verbose && println("iteration $t: Norm: $nrmB") - oldnrmB > (1-ϵ)*nrmB && break - oldnrmB = nrmB - end - for i=1:size(B.Y, 2) - scale!(view(B.Y, :, i), 1/√Λ[i]) - end - Base.LinAlg.SVD(B.X, sqrt.(Λ), B.Y') -end diff --git a/src/simple.jl b/src/simple.jl index 65b359bb..72d2f703 100644 --- a/src/simple.jl +++ b/src/simple.jl @@ -21,7 +21,7 @@ end @inline converged(p::PowerMethodIterable) = p.residual ≤ p.tol -@inline start(p::PowerMethodIterable) = 0 +@inline start(p::PowerMethodIterable) = 1 @inline done(p::PowerMethodIterable, iteration::Int) = iteration > p.maxiter || converged(p) @@ -54,12 +54,6 @@ function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxit PowerMethodIterable(A, x, tol, maxiter, zero(T), similar(x), similar(x), realmax(real(T))) end -function powm_iterable(A; kwargs...) - x0 = rand(Complex{real(eltype(A))}, size(A, 1)) - scale!(x0, one(eltype(A)) / norm(x0)) - powm_iterable!(A, x0; kwargs...) -end - """ powm(B; kwargs...) -> λ, x, [history] @@ -138,6 +132,7 @@ function powm!(B, x; for (iteration, residual) = enumerate(iterable) nextiter!(history, mvps = 1) + push!(history, :resnorm, residual) verbose && @printf("%3d\t%1.2e\n", iteration, residual) end diff --git a/src/svdl.jl b/src/svdl.jl index 040a11a7..1d948f8c 100644 --- a/src/svdl.jl +++ b/src/svdl.jl @@ -352,10 +352,10 @@ function build(log::ConvergenceHistory, A, q::AbstractVector{T}, k::Int) where { m, n = size(A) Tr = typeof(real(one(T))) β = norm(q) - @blas! q *= inv(β) + q .*= inv(β) p = A*q α = convert(Tr, norm(p)) - @blas! p *= inv(α) + p .*= inv(α) bidiag = Bidiagonal([α], Tr[], @static VERSION < v"0.7.0-DEV.884" ? true : :U) extend!(log, A, PartialFactorization(reshape(p, m, 1), reshape(q, n, 1), bidiag, β), k) end @@ -393,7 +393,7 @@ function thickrestart!(A, L::PartialFactorization{T,Tr}, #@assert ρ[i] ≈ f⋅L.P[:, i] f -= L.P*ρ α = convert(Tr, norm(f)) - @blas! f *= inv(α) + f .*= inv(α) L.P = [L.P f] g = A'f - α*L.Q[:, end] @@ -457,7 +457,7 @@ function harmonicrestart!(A, L::PartialFactorization{T,Tr}, pinv(full(L.B))*r0 else rethrow(exc) end end::Vector{Tr} - @blas! r *= L.β + r .*= L.β M::Matrix{T} = view(M,1:m, :) + r*view(M,m+1:m+1,:) M2 = zeros(T, m+1, k+1) @@ -474,16 +474,16 @@ function harmonicrestart!(A, L::PartialFactorization{T,Tr}, f = A*view(Q,:,k+1) f -= P*(P'f) α = convert(Tr, norm(f)) - @blas! f *= inv(α) + f .*= inv(α) P = [P f] B = UpperTriangular{Tr,Matrix{Tr}}([(Diagonal(Σ)*triu(R')); zeros(Tr,1,k) α]) #@assert size(P, 2) == size(B, 1) == size(Q, 2) g = A'f q = view(Q,:,k+1) #g-= (g⋅q)*q - @blas! g -= (g⋅q)*q + g .-= (g⋅q)*q β = convert(Tr, norm(g)) - @blas! g *= inv(β) + g .*= inv(β) #@assert size(P, 2) == size(Q, 2) == size(B, 2) L.P = P L.Q = Q @@ -572,7 +572,7 @@ function extend!( end β = norm(q) - @blas! q *= inv(β) + q .*= inv(β) L.Q = [L.Q q] j==k && break @@ -580,7 +580,7 @@ function extend!( log.mvps+=1 #p = A*q - β*p A_mul_B!(p, A, q) - @blas! p -= β*view(L.P, :, j) + p .-= β*view(L.P, :, j) if orthleft #Orthogonalize left Lanczos vector #Do double classical Gram-Schmidt reorthogonalization @@ -592,7 +592,7 @@ function extend!( end α = norm(p) - @blas! p *= inv(α) + p .*= inv(α) if isa(L.B, Bidiagonal) || isa(L.B, BrokenArrowBidiagonal) push!(L.B.dv, α) push!(L.B.ev, β) diff --git a/test/bicgstabl.jl b/test/bicgstabl.jl index 24168e8f..9621a05c 100644 --- a/test/bicgstabl.jl +++ b/test/bicgstabl.jl @@ -18,11 +18,18 @@ n = 20 x1, his1 = bicgstabl(A, b, l, max_mv_products = 100, log = true, tol = tol) @test isa(his1, ConvergenceHistory) @test norm(A * x1 - b) / norm(b) ≤ tol + + # With an initial guess + x_guess = rand(T, n) + x2, his2 = bicgstabl!(x_guess, A, b, l, max_mv_products = 100, log = true, tol = tol) + @test isa(his2, ConvergenceHistory) + @test x2 == x_guess + @test norm(A * x2 - b) / norm(b) ≤ tol # Do an exact LU decomp of a nearby matrix F = lufact(A + rand(T, n, n)) - x2, his2 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, tol = tol) - @test norm(A * x2 - b) / norm(b) ≤ tol + x3, his3 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, tol = tol) + @test norm(A * x3 - b) / norm(b) ≤ tol end end end diff --git a/test/chebyshev.jl b/test/chebyshev.jl index a38b3460..4168fade 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -22,20 +22,37 @@ srand(1234321) @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) A = randSPD(T, n) b = rand(T, n) - λ_min, λ_max = approx_eigenvalue_bounds(A) tol = √(eps(real(T))) - x, history = chebyshev(A, b, λ_min, λ_max, tol=tol, maxiter=10n, log=true) - @test isa(history, ConvergenceHistory) - @test history.isconverged - @test norm(A * x - b) / norm(b) ≤ tol + # Without a preconditioner + begin + λ_min, λ_max = approx_eigenvalue_bounds(A) + x0 = rand(n) + x, history = chebyshev(A, b, λ_min, λ_max, tol=tol, maxiter=10n, log=true) + @test isa(history, ConvergenceHistory) + @test history.isconverged + @test norm(A * x - b) / norm(b) ≤ tol + end - # Preconditioned solve - B = randSPD(T, n) - B_fact = cholfact!(B) - λ_min, λ_max = approx_eigenvalue_bounds(B_fact \ A) - x, history = chebyshev(A, b, λ_min, λ_max, Pl = B_fact, tol=tol, maxiter=10n, log=true) - @test history.isconverged - @test norm(A * x - b) / norm(b) ≤ tol + # With an initial guess + begin + λ_min, λ_max = approx_eigenvalue_bounds(A) + x0 = rand(T, n) + x, history = chebyshev!(x0, A, b, λ_min, λ_max, tol=tol, maxiter=10n, log=true) + @test isa(history, ConvergenceHistory) + @test history.isconverged + @test x == x0 + @test norm(A * x - b) / norm(b) ≤ tol + end + + # With a preconditioner + begin + B = randSPD(T, n) + B_fact = cholfact!(B) + λ_min, λ_max = approx_eigenvalue_bounds(B_fact \ A) + x, history = chebyshev(A, b, λ_min, λ_max, Pl = B_fact, tol=tol, maxiter=10n, log=true) + @test history.isconverged + @test norm(A * x - b) / norm(b) ≤ tol + end end end diff --git a/test/factorization.jl b/test/factorization.jl deleted file mode 100644 index d8af5908..00000000 --- a/test/factorization.jl +++ /dev/null @@ -1,11 +0,0 @@ -using IterativeSolvers -using Base.Test - -@testset "IDfact" begin - srand(1) - M = randn(4,5) - k = 3 - - F = idfact(M, k, 3) - @test vecnorm(F.B * F.P - M) ≤ 2svdvals(M)[k + 1] -end diff --git a/test/hessenberg.jl b/test/hessenberg.jl index f8fe1064..7a163db8 100644 --- a/test/hessenberg.jl +++ b/test/hessenberg.jl @@ -1,8 +1,6 @@ using IterativeSolvers using Base.Test -import IterativeSolvers: Hessenberg - @testset "Hessenberg" begin # Some well-conditioned Hessenberg matrix @@ -24,7 +22,7 @@ import IterativeSolvers: Hessenberg 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.42175+0.0im ] - for H = [H1, H2] + for H = (H1, H2) T = eltype(H) # Fist standard basis vector as rhs @@ -32,7 +30,7 @@ import IterativeSolvers: Hessenberg # Compare \ against the optimized version. solution_with_residual = copy(rhs) - A_ldiv_B!(IterativeSolvers.Hessenberg(copy(H)), solution_with_residual) + A_ldiv_B!(IterativeSolvers.FastHessenberg(copy(H)), solution_with_residual) solution = H \ rhs # First part is the solution diff --git a/test/history.jl b/test/history.jl index fa09a2af..86232133 100644 --- a/test/history.jl +++ b/test/history.jl @@ -8,6 +8,15 @@ using Base.Test RecipesBase.is_key_supported(k::Symbol) = k == :sep ? false : true + begin + history = ConvergenceHistory(partial = false) + history.iters = 3 + history.isconverged = true + @test string(history) == "Converged after 3 iterations." + history.isconverged = false + @test string(history) == "Not converged after 3 iterations." + end + # No plottables begin history = ConvergenceHistory(partial = false) diff --git a/test/lsmr.jl b/test/lsmr.jl index a0b68698..21b30a3e 100644 --- a/test/lsmr.jl +++ b/test/lsmr.jl @@ -1,7 +1,7 @@ using IterativeSolvers using Base.Test -import Base: size, A_mul_B!, Ac_mul_B!, eltype, similar, scale!, copy!, fill!, length +import Base: size, A_mul_B!, Ac_mul_B!, eltype, similar, scale!, copy!, fill!, length, broadcast! import Base.LinAlg: norm # Type used in Dampenedtest @@ -15,6 +15,12 @@ end eltype(a::DampenedVector) = promote_type(eltype(a.y), eltype(a.x)) norm(a::DampenedVector) = sqrt(norm(a.y)^2 + norm(a.x)^2) +function Base.Broadcast.broadcast!(f::Tf, to::DampenedVector, from::DampenedVector, args...) where {Tf} + to.x .= f.(from.x, args...) + to.y .= f.(from.y, args...) + to +end + function copy!(a::DampenedVector{Ty, Tx}, b::DampenedVector{Ty, Tx}) where {Ty, Tx} copy!(a.y, b.y) copy!(a.x, b.x) diff --git a/test/minres.jl b/test/minres.jl index ddad55bc..d81ee2a9 100644 --- a/test/minres.jl +++ b/test/minres.jl @@ -27,12 +27,16 @@ n = 15 @testset "Hermitian Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) A, x, b = hermitian_problem(T, n) tol = sqrt(eps(real(T))) + x0 = rand(T, n) - x_approx, hist = minres(A, b, maxiter = 10n, tol = tol, log = true) + x1, hist1 = minres(A, b, maxiter = 10n, tol = tol, log = true) + x2, hist2 = minres!(x0, A, b, maxiter = 10n, tol = tol, log = true) - @test isa(hist, ConvergenceHistory) - @test norm(b - A * x_approx) / norm(b) ≤ tol - @test hist.isconverged + @test isa(hist1, ConvergenceHistory) + @test norm(b - A * x1) / norm(b) ≤ tol + @test hist1.isconverged + @test norm(b - A * x2) / norm(b) ≤ tol + @test x2 == x0 end @testset "Skew-Hermitian Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128) diff --git a/test/rlinalg.jl b/test/rlinalg.jl deleted file mode 100644 index b730d00f..00000000 --- a/test/rlinalg.jl +++ /dev/null @@ -1,31 +0,0 @@ -using IterativeSolvers -using Base.Test - -@testset "Randomized Linear Algebra" begin - -srand(1234321) - -m = n = 100 -B = randn(m, n) -nB = norm(B) -p = 1e-5 #Probability of failure - -for j = 1 : n - @test rnorm(B, 2j, p) ≥ nB - @test rnorms(B, j, p) ≥ nB -end - -A = B * B' -k = 1 -l, u = reigmin(A, k, p) -@test l ≤ eigmin(A) ≤ u - -l, u = reigmax(A, k, p) -@test l ≤ eigmax(A) ≤ u - -l, u = rcond(A, k, p) -@test l ≤ cond(A) ≤ u - -@test_throws ArgumentError IterativeSolvers.randnn(Char, m) -@test_throws ArgumentError IterativeSolvers.randnn(Char, m, n) -end diff --git a/test/rsvd.jl b/test/rsvd.jl deleted file mode 100644 index ddfc696b..00000000 --- a/test/rsvd.jl +++ /dev/null @@ -1,60 +0,0 @@ -using IterativeSolvers -using Base.Test - -@testset "rsvd" begin - -srand(1234321) - -@testset "Small wide rectangular" begin - A = [1. 2 3 4; 5 6 7 8] - S1 = svdfact(A) - S2 = rsvdfact(A, 2, 0) - - @test vecnorm(abs.(S1[:U]) - abs.(S2[:U])) ≤ √(eps()) - @test vecnorm(abs.(S1[:Vt]) - abs.(S2[:Vt])) ≤ √(eps()) - @test norm(S1[:S] - S2[:S]) ≤ √(eps()) -end - -@testset "Small tall rectangular" begin - A = [1. 2; 3 4; 5 6; 7 8] - S1 = svdfact(A) - S2 = rsvdfact(A, 2, 0) - - @test vecnorm(abs.(S1[:U]) - abs.(S2[:U])) ≤ √(eps()) - @test vecnorm(abs.(S1[:Vt]) - abs.(S2[:Vt])) ≤ √(eps()) - @test norm(S1[:S] - S2[:S]) ≤ √(eps()) -end - -@testset "Small square" begin - A = [1. 2 3; 4 5 6; 7 8 9] - S1 = svdfact(A) - S2 = rsvdfact(A, 3, 0) - - @test vecnorm(abs.(S1[:U]) - abs.(S2[:U])) ≤ √(eps()) - @test vecnorm(abs.(S1[:Vt]) - abs.(S2[:Vt])) ≤ √(eps()) - @test norm(S1[:S] - S2[:S]) ≤ √(eps()) -end - -@testset "Low rank" begin #Issue #42 - n = 10 - r = 2 - A = randn(n, r) * randn(r, n) - S = svdvals(A) - for nvals = 1 : r - S1 = IterativeSolvers.rsvdvals(A, nvals, r - nvals) - for i = 1 : nvals - @test abs(S[i] - S1[i]) ≤ n^2 * r * eps() - end - end -end - -@testset "rrange_adaptive" begin - A = [1. 2 3; 4 5 6; 7 8 9] - @test size(IterativeSolvers.rrange_adaptive(A, 3, 1e-3)) == (3,2) -end - -@testset "rrange" begin - A = [1. 2 3; 4 5 6; 7 8 9] - @test_throws ArgumentError IterativeSolvers.rrange(A, 20) -end -end diff --git a/test/rsvd_fnkz.jl b/test/rsvd_fnkz.jl deleted file mode 100644 index f79fd049..00000000 --- a/test/rsvd_fnkz.jl +++ /dev/null @@ -1,28 +0,0 @@ -using IterativeSolvers -using Base.Test - -@testset "rsvd_fnkz" begin - -srand(1234321) - -@testset "Rank 2 matrix" begin - A = reshape(1 : 64, 8, 8) - S = rsvd_fnkz(A, 2, ϵ=1e-9) - @test vecnorm(A - S[:U] * Diagonal(S[:S]) * S[:Vt]) ≤ 1e-9 -end - -@testset "Linearly dependent matrix" begin - A = [collect(1 : 8) collect(1 : 8) collect(1 : 8) collect(1 : 8)] - S = rsvd_fnkz(A, 4) - @test vecnorm(A - S[:U] * Diagonal(S[:S]) * S[:Vt]) ≤ 1e-9 -end - -#SVD of a matrix of zeros should be empty -@testset "Zero matrix" begin - S = rsvd_fnkz(zeros(10, 10), 1) - @test S[:U] == zeros(10, 0) - @test S[:S] == [] - @test S[:Vt] == zeros(0, 10) -end - -end diff --git a/test/runtests.jl b/test/runtests.jl index 560130ba..056880f2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,11 +37,6 @@ include("svdl.jl") include("lsqr.jl") include("lsmr.jl") -#Randomized algorithms -include("rlinalg.jl") -include("rsvd.jl") -include("rsvd_fnkz.jl") - #History data structure include("history.jl")