Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

An implementation of DIOM #63

Merged
merged 9 commits into from
Oct 23, 2018
Prev Previous commit
Next Next commit
Use exact residual norm in DIOM
  • Loading branch information
amontoison committed Oct 14, 2018
commit 1e5569b346c4580c96da185768bf640d6eebc172
73 changes: 39 additions & 34 deletions src/diom.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,16 +54,16 @@ function diom{T <: Number}(A :: AbstractLinearOperator, b :: AbstractVector{T};

# Set up workspace.
mem = min(memory, itmax) # Memory.
P = zeros(n, mem) # Directions for x.
V = zeros(n, mem) # Preconditioned Krylov vectors.
V = zeros(n, mem) # Preconditioned Krylov vectors, orthogonal basis for {r₀, AM⁻¹r₀, (AM⁻¹)²r₀, ..., (AM⁻¹)ᵐ⁻¹r₀}.
P = zeros(n, mem) # Directions for x : Pₘ = Vₘ(Uₘ)⁻¹.
H = zeros(mem+2) # Last column of the band hessenberg matrix Hₘ = LₘUₘ.
# Each column has at most mem + 2 nonzero elements. hᵢ.ₘ is stored as H[m-i+2].
# Each column has at most mem + 1 nonzero elements. hᵢ.ₘ is stored as H[m-i+2].
# m-i+2 represents the indice of the diagonal where hᵢ.ₘ is located.
# In addition of that, the last column of Uₘ is stored in H.
L = zeros(mem) # Last mem Pivots of Lₘ.
p = zeros(Bool, mem) # Last mem permutations.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe p could be a BitArray?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I didn't know that a Bool was an Uint8 in Julia and not a single bit.


# Initial ξ and V.
# Initial ξ and V.
ξ = rNorm
V[:,1] = r / rNorm

Expand All @@ -77,61 +77,56 @@ function diom{T <: Number}(A :: AbstractLinearOperator, b :: AbstractVector{T};
# Update iteration index.
iter = iter + 1

# Set position in circular stack where iter-th Krylov vector should go.
pos = mod(iter-1, mem) + 1 # Position corresponding to iter in the circular stack.
next_pos = mod(iter, mem) + 1 # Position corresponding to iter + 1 in the circular stack.
# Set position in circulars stacks.
pos = mod(iter-1, mem) + 1 # Position corresponding to pₘ and vₘ in circular stacks P and V.
next_pos = mod(iter, mem) + 1 # Position corresponding to vₘ₊₁ in the circular stack V.

# Incomplete Arnoldi procedure.
z = M * V[:,pos] # Forms pₘ
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder why Github isn't rendering UTF-8 characters correctly. Is your editor set to use UTF-8 encoding?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, my editor is set to use UTF-8 encoding. Github is rendering UTF-8 correctly for me, maybe it's your web browser?
encodage

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both Safari and Chrome seem to be missing characters...

w = A * z # Forms vₘ₊₁
# hₘ₋ₘₑₘ.ₘ = 0
H[mem+2] = 0
for i = max(1, iter-mem+1) : iter
ipos = mod(i-1, mem) + 1 # Position of vᵢ in the circular stack
jpos = iter - i + 2
H[jpos] = @kdot(n, w, V[:,ipos]) # hᵢ.ₘ = < A * vₘ , vᵢ >
@kaxpy!(n, -H[jpos], V[:,ipos], w) # w ← w - hᵢ.ₘ * vᵢ
ipos = mod(i-1, mem) + 1 # Position corresponding to vᵢ in the circular stack V.
diag = iter - i + 2
H[diag] = @kdot(n, w, V[:,ipos]) # hᵢ.ₘ = < A * vₘ , vᵢ >
@kaxpy!(n, -H[diag], V[:,ipos], w) # w ← w - hᵢ.ₘ * vᵢ
end
# Compute hₘ₊₁.ₘ and vₘ₊₁.
H[1] = @knrm2(n, w) # hₘ₊₁.ₘ = ‖vₘ₊₁‖
if H[1] ≉ 0 # hₘ₊₁.ₘ ≈ 0 ⇒ "lucky breakdown"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't there a tolerance involved in ? In my experience, we want to test for != 0.

Copy link
Member Author

@amontoison amontoison Oct 22, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, there is a tolerance involved in , x ≈ y return true if norm(x-y) ≤ rtol * max(|x|,|y|). With rtol that depend of the type of x and y. It's sqrt(eps) with eps the precision machine of Float16, Float32, Float64....

But when we do x ≈ 0, It's mean that |x| ≤ rtol * |x| ↔ x = 0. I used this notation because I find it nicer than ==. Also I like this notation because it means that x is a zero for the computer in the working precision but it's possible that it's not a zero in another precision or in exact arithmetic.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK.

V[:,next_pos] = w / H[1] # vₘ₊₁ = w / hₘ₊₁.ₘ
end
# It's possible that uₘ₋ₘₑₘ.ₘ ≠ 0 when m ≥ mem + 1
if iter ≥ mem + 2
H[mem+2] = 0 # hₘ₋ₘₑₘ.ₘ = 0
end

# Compute LU factorization with partial pivoting of Hₘ by computing the last column of Uₘ.
# Update the LU factorization with partial pivoting of H.
# Compute the last column of Uₘ.
if iter ≥ 2
for i = max(2,iter-mem+1) : iter
lpos = mod(i-1, mem) + 1 # Position of lᵢ.ᵢ₋₁ in the circular stack
jpos = iter - i + 2
lpos = mod(i-1, mem) + 1 # Position corresponding to lᵢ.ᵢ₋₁ in the circular stack L.
diag = iter - i + 2
next_diag = diag + 1
if p[lpos]
# row i-1 and i are permuted
H[jpos], H[jpos+1] = H[jpos+1], H[jpos]
# The rows i-1 and i are permuted.
H[diag], H[next_diag] = H[next_diag], H[diag]
end
# uᵢ.ₘ ← hᵢ.ₘ - lᵢ.ᵢ₋₁ * uᵢ₋₁.ₘ
H[jpos] = H[jpos] - L[lpos] * H[jpos+1]
H[diag] = H[diag] - L[lpos] * H[next_diag]
end
# Compute ξₘ the last composant of zₘ = β(Lₘ)⁻¹e₁.
# Compute ξₘ the last component of zₘ = β(Lₘ)⁻¹e₁.
if !p[pos] # p[pos] ⇒ ξₘ = ξₘ₋₁
# ξₘ = -lₘ.ₘ₋₁ * ξₘ₋₁
ξ = - L[pos] * ξ
end
end

# Update residual norm estimate.
# ‖ b - Axₘ ‖ = hₘ₊₁.ₘ * |ξₘ / uₘ.ₘ|
rNorm = H[1] * abs(ξ / H[2])
push!(rNorms, rNorm)

# Update stopping criterion.
solved = rNorm ≤ ε
tired = iter ≥ itmax

# Compute the direction pₘ, the last column of Pₘ = Vₘ(Uₘ)⁻¹.
for i = max(1,iter-mem) : iter - 1
ipos = mod(i-1, mem) + 1
jpos = iter - i + 2
ipos = mod(i-1, mem) + 1 # Position corresponding to pᵢ in the circular stack P.
diag = iter - i + 2
# z ← z - uᵢ.ₘ * pᵢ
@kaxpy!(n, -H[jpos], P[:,ipos], z)
@kaxpy!(n, -H[diag], P[:,ipos], z)
end

# Determine if interchange between hₘ₊₁.ₘ and uₘ.ₘ is needed and compute next pivot lₘ₊₁.ₘ.
Expand All @@ -152,17 +147,27 @@ function diom{T <: Number}(A :: AbstractLinearOperator, b :: AbstractVector{T};
# Compute solution xₘ.
if p[pos]
# xₘ = xₘ₋ₙ + ξₘ₋ₙ * pₘ
# m-n is the last iteration without permutation at the next step
# x_old = xₘ₋ₙ, with m-n is the last iteration without permutation at the next step
x = x_old + ξ * P[:,pos]
else
# xₘ = xₘ₋₁ + ξₘ * pₘ
@kaxpy!(n, ξ, P[:,pos], x)
end

# Update x_old.
# Update x_old and residual norm.
if !p[next_pos]
x_old = copy(x)
# ‖ b - Axₘ ‖ = hₘ₊₁.ₘ * |ξₘ / uₘ.ₘ| without pivoting
rNorm = H[1] * abs(ξ / H[2])
else
# ‖ b - Axₘ ‖ = hₘ₊₁.ₘ * |ξₘ / hₘ₊₁.ₘ| = |ξₘ| with pivoting
rNorm = abs(ξ)
end
push!(rNorms, rNorm)

# Update stopping criterion.
solved = rNorm ≤ ε
tired = iter ≥ itmax
verbose && @printf("%5d %7.1e\n", iter, rNorm)
end
verbose && @printf("\n")
Expand Down