diff --git a/src/Krylov.jl b/src/Krylov.jl
index da80eb7f0..8268144d4 100644
--- a/src/Krylov.jl
+++ b/src/Krylov.jl
@@ -77,6 +77,7 @@ include("cg.jl")
 include("cg_lanczos.jl")
 include("minres.jl")
 include("dqgmres.jl")
+include("diom.jl")
 include("symmlq.jl")
 include("cr.jl")
 
diff --git a/src/diom.jl b/src/diom.jl
new file mode 100644
index 000000000..b90095375
--- /dev/null
+++ b/src/diom.jl
@@ -0,0 +1,173 @@
+# An implementation of DIOM for the solution of the square linear system Ax = b.
+#
+# This method is described in
+#
+# Y. Saad, Iterative methods for sparse linear systems.
+# PWS Publishing Company, Boston, USA, 1996.
+#
+# Y. Saad, Practical use of some krylov subspace methods for solving indefinite and nonsymmetric linear systems.
+# SIAM journal on scientific and statistical computing, 5(1), pp. 203--228, 1984.
+# 
+# Alexis Montoison, <alexis.montoison@polymtl.ca>
+# Montreal, September 2018.
+
+export diom
+
+"""Solve the consistent linear system Ax = b using direct incomplete orthogonalization method.
+
+DIOM is similar to CG with partial reorthogonalization.
+
+An advantage of DIOM is that nonsymmetric or symmetric indefinite or both nonsymmetric 
+and indefinite systems of linear equations can be handled by this single algorithm.
+
+This implementation allows a right preconditioning with M.
+"""
+function diom{T <: Number}(A :: AbstractLinearOperator, b :: AbstractVector{T};
+                              M :: AbstractLinearOperator=opEye(size(A,1)),
+                              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
+                              itmax :: Int=0, memory :: Int=20, verbose :: Bool=false)
+
+  m, n = size(A)
+  m == n || error("System must be square")
+  size(b, 1) == m || error("Inconsistent problem size")
+  verbose && @printf("DIOM: system of size %d\n", n)
+
+  # Initial solution x₀ and residual r₀.
+  x = zeros(T, n)
+  x_old = copy(x)
+  r = copy(b)
+  # Compute β.
+  rNorm = @knrm2(n, r) # rNorm = ‖r₀‖
+  rNorm ≈ 0 && return x, SimpleStats(true, false, [rNorm], [], "x = 0 is a zero-residual solution")
+
+  iter = 0
+  itmax == 0 && (itmax = 2*n)
+
+  rNorms = [rNorm;]
+  ε = atol + rtol * rNorm
+  verbose && @printf("%5d  %7.1e\n", iter, rNorm)
+
+  # Set up workspace.
+  mem = min(memory, itmax) # Memory.
+  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 + 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 = BitArray(mem) # Last mem permutations.
+
+  # Initial ξ₁ and V₁.
+  ξ = rNorm
+  V[:,1] = r / rNorm
+
+  # Stopping criterion.
+  solved = rNorm ≤ ε
+  tired = iter ≥ itmax
+  status = "unknown"
+
+  while !(solved || tired)
+
+    # Update iteration index.
+    iter = iter + 1
+
+    # 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ₘ
+    w = A * z # Forms vₘ₊₁
+    for i = max(1, iter-mem+1) : iter
+      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"
+      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
+
+    # 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 corresponding to lᵢ.ᵢ₋₁ in the circular stack L.
+        diag = iter - i + 2
+        next_diag = diag + 1
+        if p[lpos]
+          # The rows i-1 and i are permuted.
+          H[diag], H[next_diag] = H[next_diag], H[diag]
+        end
+        # uᵢ.ₘ ← hᵢ.ₘ - lᵢ.ᵢ₋₁ * uᵢ₋₁.ₘ
+        H[diag] = H[diag] - L[lpos] * H[next_diag]
+      end
+      # Compute ξₘ the last component of zₘ = β(Lₘ)⁻¹e₁.
+      if !p[pos] # p[pos] ⇒ ξₘ = ξₘ₋₁
+        # ξₘ = -lₘ.ₘ₋₁ * ξₘ₋₁
+        ξ = - L[pos] * ξ
+      end
+    end
+
+    # 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 # Position corresponding to pᵢ in the circular stack P.
+      diag = iter - i + 2
+      # z ← z - uᵢ.ₘ * pᵢ
+      @kaxpy!(n, -H[diag], P[:,ipos], z)
+    end
+
+    # Determine if interchange between hₘ₊₁.ₘ and uₘ.ₘ is needed and compute next pivot lₘ₊₁.ₘ.
+    if abs(H[2]) < H[1]
+      p[next_pos] = true
+      # pₘ = z / hₘ₊₁.ₘ
+      P[:,pos] = z / H[1]
+      # lₘ₊₁.ₘ = uₘ.ₘ / hₘ₊₁.ₘ
+      L[next_pos] = H[2] / H[1]
+    else
+      p[next_pos] = false
+      # pₘ = z / uₘ.ₘ
+      P[:,pos] = z / H[2]
+      # lₘ₊₁.ₘ = hₘ₊₁.ₘ / uₘ.ₘ
+      L[next_pos] = H[1] / H[2]
+    end
+
+    # Compute solution xₘ.
+    if p[pos]
+      # xₘ = xₘ₋ₙ + ξₘ₋ₙ * pₘ
+      # 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 and residual norm.
+    if !p[next_pos]
+      copy!(x_old, 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")
+
+  status = tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol"
+  stats = SimpleStats(solved, false, rNorms, T[], status)
+  return (x, stats)
+end
diff --git a/src/variants.jl b/src/variants.jl
index 1b349662a..c1c3abdbb 100644
--- a/src/variants.jl
+++ b/src/variants.jl
@@ -7,7 +7,7 @@ function preallocated_LinearOperator(A)
 end
 
 # Variants where A is a matrix without specific properties
-for fn in (:cgls, :cgne, :craig, :craigmr, :crls, :crmr, :lslq, :lsmr, :lsqr, :dqgmres)
+for fn in (:cgls, :cgne, :craig, :craigmr, :crls, :crmr, :lslq, :lsmr, :lsqr, :dqgmres, :diom)
   @eval begin
     # Variant for A given as a dense array and b given as a dense vector
     $fn{TA <: Number, Tb <: Number}(A :: Array{TA,2}, b :: Vector{Tb}, args...; kwargs...) =
diff --git a/test/runtests.jl b/test/runtests.jl
index 20aa3d0f2..2660f7199 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,7 @@ using Base.Test
 using Krylov
 using LinearOperators
 
+include("test_diom.jl")
 include("test_dqgmres.jl")
 include("gen_lsq.jl")
 include("check_min_norm.jl")
diff --git a/test/test_diom.jl b/test/test_diom.jl
new file mode 100644
index 000000000..c21f26070
--- /dev/null
+++ b/test/test_diom.jl
@@ -0,0 +1,89 @@
+include("get_div_grad.jl")
+
+diom_tol = 1.0e-6
+
+# Symmetric and positive definite systems (cubic spline matrix).
+n = 10
+A = spdiagm((ones(n-1), 4*ones(n), ones(n-1)), (-1, 0, 1))
+b = A * [1:n;]
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Symmetric indefinite variant.
+A = A - 3 * speye(n)
+b = A * [1:n;]
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Nonsymmetric and positive definite systems.
+A = [i == j ? 10.0 : i < j ? 1.0 : -1.0 for i=1:n, j=1:n]
+b = A * [1:n;]
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Nonsymmetric indefinite variant.
+A = [i == j ? 10*(-1.0)^(i*j) : i < j ? 1.0 : -1.0 for i=1:n, j=1:n]
+b = A * [1:n;]
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Code coverage.
+(x, stats) = diom(full(A), b)
+show(stats)
+
+# Sparse Laplacian.
+A = get_div_grad(16, 16, 16)
+b = ones(size(A, 1))
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Symmetric indefinite variant, almost singular.
+A = A - 5 * speye(size(A, 1))
+(x, stats) = diom(A, b)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
+
+# Test b == 0
+(x, stats) = diom(A, zeros(size(A,1)))
+@test x == zeros(size(A,1))
+@test stats.status == "x = 0 is a zero-residual solution"
+
+# Test integer values
+A = spdiagm((ones(Int, n-1), 4*ones(Int, n), ones(Int, n-1)), (-1, 0, 1)); b = A * [1:n;]
+(x, stats) = diom(A, b)
+@test stats.solved
+
+# Test with Jacobi (or diagonal) preconditioner
+A = ones(10,10) + 9 * eye(10)
+b = 10 * [1:10;]
+M = 1/10 * opEye(10)
+(x, stats) = diom(A, b, M=M)
+show(stats)
+r = b - A * x
+resid = norm(r) / norm(b)
+@printf("DIOM: Relative residual: %8.1e\n", resid)
+@test(resid ≤ diom_tol)
+@test(stats.solved)
diff --git a/test/test_variants.jl b/test/test_variants.jl
index 7c79e502d..9bab5f2b5 100644
--- a/test/test_variants.jl
+++ b/test/test_variants.jl
@@ -1,5 +1,5 @@
 # Tests of variants.jl
-for fn in (:cg_lanczos, :cg_lanczos_shift_seq, :cg, :cgls, :cgne, :cr, :craig, :craigmr, :crls, :crmr, :lslq, :lsmr, :lsqr, :minres, :symmlq, :dqgmres)
+for fn in (:cg_lanczos, :cg_lanczos_shift_seq, :cg, :cgls, :cgne, :cr, :craig, :craigmr, :crls, :crmr, :lslq, :lsmr, :lsqr, :minres, :symmlq, :dqgmres, :diom)
   for TA in (Int32, Int64, Float32, Float64)
     for IA in (Int32, Int64)
       for Tb in (Int32, Int64, Float32, Float64)