Skip to content

Commit

Permalink
added residual smoothing as in the python code https://github.com/ast…
Browse files Browse the repository at this point in the history 
…udillor/idrs_py/blob/master/idrs.py
  • Loading branch information
@jiahao
RomainBrault authored and jiahao committed Aug 23, 2016
1 parent e38feea commit 6f4df86
Show file tree
Hide file tree
Showing 3 changed files with 111 additions and 59 deletions.
108 changes: 72 additions & 36 deletions src/idrs.jl
View file
@@ -1,3 +1,4 @@
import Base.LinAlg.BLAS: axpy!
export idrs, idrs!

####
Expand All @@ -24,9 +25,9 @@ IDRsSolver.jl
The Induced Dimension Reduction method is a family of simple and fast Krylov
subspace algorithms for solving large nonsymmetric linear systems. The idea
behind the IDR(s) variant is to generate residuals that are in the nested
subspaces of shrinking dimension s.
subspaces of shrinking dimension s.
Syntax
------
Expand All @@ -39,83 +40,91 @@ Syntax
Arguments
---------
s -- dimension reduction number. Normally, a higher s gives faster convergence,
s -- dimension reduction number. Normally, a higher s gives faster convergence,
but also makes the method more expensive.
tol -- tolerance of the method.
tol -- tolerance of the method.
maxiter -- maximum number of iterations
smoothing -- (Logical) specifies if residual smoothing must be applied. Default is false.
Output
------
X -- Approximated solution by IDR(s)
h -- Convergence history
The [`ConvergenceHistory`](https://github.com/JuliaLang/IterativeSolvers.jl/issues/6) type provides information about the iteration history.
The [`ConvergenceHistory`](https://github.com/JuliaLang/IterativeSolvers.jl/issues/6) type provides information about the iteration history.
- `isconverged::Bool`, a flag for whether or not the algorithm is converged.
- `threshold`, the convergence threshold
- `residuals::Vector`, the value of the convergence criteria at each iteration
- `residuals::Vector`, the value of the convergence criteria at each iteration
References
----------
[1] IDR(s): a family of simple and fast algorithms for solving large
[1] IDR(s): a family of simple and fast algorithms for solving large
nonsymmetric linear systems. P. Sonneveld and M. B. van Gijzen
SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035--1062, 2008
[2] Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits
SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035--1062, 2008
[2] Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits
Bi-orthogonality Properties. M. B. van Gijzen and P. Sonneveld
ACM Trans. Math. Software,, Vol. 38, No. 1, pp. 5:1-5:19, 2011
[3] This file is a translation of the following MATLAB implementation:
http://ta.twi.tudelft.nl/nw/users/gijzen/idrs.m
[4] IDR(s)' webpage http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html
"""
idrs(A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter = length(b)^2) =
idrs_core!(zerox(A,b), linsys_op, (A,), b, s, tol, maxiter)
idrs(A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter = length(b)^2, smoothing=false) =
idrs_core!(zerox(A, b), linsys_op, (A,), b, s, tol, maxiter; smoothing=smoothing)

idrs!(x, A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter=length(x)^2) =
idrs_core!(x, linsys_op, (A,), b, s, tol, maxiter)
idrs!(x, A, b; s = 8, tol=sqrt(eps(typeof(real(b[1])))), maxiter=length(x)^2, smoothing=false) =
idrs_core!(x, linsys_op, (A,), b, s, tol, maxiter; smoothing=smoothing)

function idrs_core!{T}(X, op, args, C::T,
s::Number, tol::Number, maxiter::Number; smoothing::Bool=false)

function idrs_core!{T}(X, op, args, C::T, s, tol, maxiter)
R = C - op(X, args...)::T
normR = vecnorm(R)
res = typeof(tol)[normR]
iter = 0
iter = 0

if smoothing
X_s = copy(X)
R_s = copy(R)
T_s = zeros(R)
end

if normR <= tol # Initial guess is a good enough solution
return X, ConvergenceHistory(0<= res[end] < tol, tol, length(res), res)
end

Z = zero(C)

P = T[rand!(copy(C)) for k in 1:s]
U = T[copy(Z) for k in 1:s]
G = T[copy(Z) for k in 1:s]
Q = copy(Z)
V = copy(Z)

M = eye(eltype(C),s,s)
f = zeros(eltype(C),s)
c = zeros(eltype(C),s)

om::eltype(C) = 1
iter = 0
while normR > tol && iter < maxiter
for i in 1:s,
f[i] = vecdot(P[i], R)
end
for k in 1:s
for k in 1:s

# Solve small system and make v orthogonal to P

c = LowerTriangular(M[k:s,k:s])\f[k:s]
copy!(V, G[k])
scale!(c[1], V)

copy!(Q, U[k])
scale!(c[1], Q)
for i = k+1:s
Expand All @@ -128,10 +137,10 @@ function idrs_core!{T}(X, op, args, C::T, s, tol, maxiter)
#V = R - V
scale!(-1., V)
axpy!(1., R, V)

copy!(U[k], Q)
axpy!(om, V, U[k])
G[k] = op(U[k], args...)
G[k] = op(U[k], args...)

# Bi-orthogonalise the new basis vectors

Expand All @@ -144,41 +153,68 @@ function idrs_core!{T}(X, op, args, C::T, s, tol, maxiter)
# New column of M = P'*G (first k-1 entries are zero)

for i in k:s
M[i,k] = vecdot(P[i],G[k])
M[i,k] = vecdot(P[i],G[k])
end

# Make r orthogonal to q_i, i = 1..k
# Make r orthogonal to q_i, i = 1..k

beta = f[k]/M[k,k]
axpy!(-beta, G[k], R)
axpy!(beta, U[k], X)

normR = vecnorm(R)
if smoothing
# T_s = R_s - R
copy!(T_s, R_s); axpy!(-1., R, T_s)

gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s)

# R_s = R_s - gamma*T_s
axpy!(-gamma, T_s, R_s)
# X_s = X_s - gamma*(X_s - X)
copy!(T_s, X_s); axpy!(-1., X, T_s); axpy!(-gamma, T_s, X_s)

normR = vecnorm(R_s)
end
res = [res; normR]
iter += 1
if normR < tol || iter > maxiter
return X, ConvergenceHistory(normR < tol, tol, length(res), res)
end
if k < s
end
if k < s
f[k+1:s] = f[k+1:s] - beta*M[k+1:s,k]
end

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

copy!(V, R)
Q = op(V, args...)::T
om = omega(Q, R)
axpy!(-om, Q, R)
axpy!(om, V, X)


normR = vecnorm(R)
if smoothing
# T_s = R_s - R
copy!(T_s, R_s); axpy!(-1., R, T_s)

gamma = vecdot(R_s, T_s)/vecdot(T_s, T_s)

# R_s = R_s - gamma*T_s
axpy!(-gamma, T_s, R_s)
# X_s = X_s - gamma*(X_s - X)
copy!(T_s, X_s); axpy!(-1., X, T_s); axpy!(-gamma, T_s, X_s)

normR = vecnorm(R_s)
end
iter += 1
res = [res; normR]
end
if smoothing
copy!(X, X_s)
end
return X, ConvergenceHistory(res[end]<tol, tol, length(res), res)
end

44 changes: 22 additions & 22 deletions src/lsmr.jl
View file
@@ -1,6 +1,6 @@
export lsmr, lsmr!

using Base.LinAlg
using Base.LinAlg

##############################################################################
## LSMR
Expand All @@ -10,7 +10,7 @@ using Base.LinAlg
## Adapted from the BSD-licensed Matlab implementation at
## http://web.stanford.edu/group/SOL/software/lsmr/
##
## A is a StridedVecOrMat or anything that implements
## A is a StridedVecOrMat or anything that implements
## A_mul_B!(α, A, b, β, c) updates c -> α Ab + βc
## Ac_mul_B!(α, A, b, β, c) updates c -> α A'b + βc
## eltype(A)
Expand Down Expand Up @@ -39,8 +39,8 @@ using Base.LinAlg
## x is initial x0. Transformed in place to the solution.
## b equals initial b. Transformed in place
## v, h, hbar are storage arrays of length size(A, 2)
function lsmr!(x, A, b, v, h, hbar;
atol::Number = 1e-6, btol::Number = 1e-6, conlim::Number = 1e8,
function lsmr!(x, A, b, v, h, hbar;
atol::Number = 1e-6, btol::Number = 1e-6, conlim::Number = 1e8,
maxiter::Integer = max(size(A,1), size(A,2)), λ::Number = 0)

# Sanity-checking
Expand Down Expand Up @@ -93,13 +93,13 @@ function lsmr!(x, A, b, v, h, hbar;

# Items for use in stopping rules.
normb = β
istop = 0
istop = 0
normr = β
normAr = α * β
tests = Tuple{Tr, Tr, Tr}[]
iter = 0
# Exit if b = 0 or A'b = 0.
if normAr != 0
if normAr != 0
while iter < maxiter
iter += 1
A_mul_B!(1, A, v, -α, u)
Expand All @@ -110,20 +110,20 @@ function lsmr!(x, A, b, v, h, hbar;
α = norm(v)
α > 0 && scale!(v, inv(α))
end

# Construct rotation Qhat_{k,2k+1}.
αhat = hypot(αbar, λ)
chat = αbar / αhat
shat = λ / αhat

# Use a plane rotation (Q_i) to turn B_i to R_i.
ρold = ρ
ρ = hypot(αhat, β)
c = αhat / ρ
s = β / ρ
θnew = s * α
αbar = c * α

# Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar.
ρbarold = ρbar
ζold = ζ
Expand All @@ -134,28 +134,28 @@ function lsmr!(x, A, b, v, h, hbar;
sbar = θnew / ρbar
ζ = cbar * ζbar
ζbar = - sbar * ζbar

# Update h, h_hat, x.
scale!(hbar, - θbar * ρ / (ρold * ρbarold))
axpy!(1, h, hbar)
axpy!(ζ / (ρ * ρbar), hbar, x)
scale!(h, - θnew / ρ)
axpy!(1, v, h)

##############################################################################
##
## Estimate of ||r||
##
##############################################################################

# Apply rotation Qhat_{k,2k+1}.
βacute = chat * βdd
βcheck = - shat * βdd

# Apply rotation Q_{k,k+1}.
βhat = c * βacute
βdd = - s * βacute

# Apply rotation Qtilde_{k-1}.
θtildeold = θtilde
ρtildeold = hypot(ρdold, θbar)
Expand All @@ -164,34 +164,34 @@ function lsmr!(x, A, b, v, h, hbar;
θtilde = stildeold * ρbar
ρdold = ctildeold * ρbar
βd = - stildeold * βd + ctildeold * βhat

τtildeold = (ζold - θtildeold * τtildeold) / ρtildeold
τd = (ζ - θtilde * τtildeold) / ρdold
d += abs2(βcheck)
normr = sqrt(d + abs2(βd - τd) + abs2(βdd))

# Estimate ||A||.
normA2 += abs2(β)
normA = sqrt(normA2)
normA2 += abs2(α)

# Estimate cond(A).
maxrbar = max(maxrbar, ρbarold)
if iter > 1
if iter > 1
minrbar = min(minrbar, ρbarold)
end
condA = max(maxrbar, ρtemp) / min(minrbar, ρtemp)

##############################################################################
##
## Test for convergence
##
##############################################################################

# Compute norms for convergence testing.
normAr = abs(ζbar)
normx = norm(x)

# Now use these norms to estimate certain other quantities,
# some of which will be small near a solution.
test1 = normr / normb
Expand All @@ -200,7 +200,7 @@ function lsmr!(x, A, b, v, h, hbar;
push!(tests, (test1, test2, test3))

t1 = test1 / (one(Tr) + normA * normx / normb)
rtol = btol + atol * normA * normx / normb
rtol = btol + atol * normA * normx / normb
# The following tests guard against extremely small values of
# atol, btol or ctol. (The user may have set any or all of
# the parameters atol, btol, conlim to 0.)
Expand Down

0 comments on commit 6f4df86

Please sign in to comment.