diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index db85283..3713c1b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,9 +15,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
-          - 'nightly'
         os:
           - ubuntu-latest
         arch:
diff --git a/Project.toml b/Project.toml
index d6e841a..aa9a727 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,17 +1,23 @@
 name = "WoodburyMatrices"
 uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
-version = "0.5.5"
+version = "0.5.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-julia = "1"
+Aqua = "0.8"
+Random = "1"
+LinearAlgebra = "1"
+SparseArrays = "1"
+Test = "1"
+julia = "1.6"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Random", "Test"]
+test = ["Aqua", "Random", "Test"]
diff --git a/src/WoodburyMatrices.jl b/src/WoodburyMatrices.jl
index 29383c4..3ce578b 100644
--- a/src/WoodburyMatrices.jl
+++ b/src/WoodburyMatrices.jl
@@ -12,13 +12,6 @@ abstract type AbstractWoodbury{T} <: Factorization{T} end
 safeinv(A) = inv(A)
 safeinv(A::SparseMatrixCSC) = safeinv(Matrix(A))
 
-myldiv!(A, B)       = ldiv!(A, B)
-myldiv!(dest, A, B) = ldiv!(dest, A, B)
-if VERSION <= v"1.4.0-DEV.635"
-    myldiv!(A::Diagonal, B)       = (B .= A.diag .\ B)
-    myldiv!(dest, A::Diagonal, B) = (dest .= A.diag .\ B)
-end
-
 include("woodbury.jl")
 include("symwoodbury.jl")
 include("sparsefactors.jl")
@@ -64,6 +57,7 @@ function _ldiv(W::AbstractWoodbury, R::AbstractMatrix)
 end
 
 \(W::AbstractWoodbury, R::AbstractMatrix) = _ldiv(W, R)
+\(W::AbstractWoodbury{T}, R::Matrix{Complex{T}}) where T<:Union{Float32, Float64} = _ldiv(W, R)  # ambiguity resolution
 \(W::AbstractWoodbury, D::Diagonal) = _ldiv(W, D)
 
 ldiv!(W::AbstractWoodbury, B::AbstractVector) = ldiv!(B, W, B)
@@ -76,18 +70,20 @@ function ldiv!(dest::AbstractVector, W::AbstractWoodbury, B::AbstractVector)
     return dest
 end
 
-function _ldiv!(dest, W, A::Union{Factorization,Diagonal}, B)
-    myldiv!(W.tmpN1, A, B)
+function _ldiv!(dest, W::AbstractWoodbury, A::Union{Factorization,Diagonal}, B)
+    ldiv!(W.tmpN1, A, B)
     mul!(W.tmpk1, W.V, W.tmpN1)
     mul!(W.tmpk2, W.Cp, W.tmpk1)
     mul!(W.tmpN2, W.U, W.tmpk2)
-    myldiv!(A, W.tmpN2)
+    ldiv!(A, W.tmpN2)
     for i = 1:length(W.tmpN2)
         @inbounds dest[i] = W.tmpN1[i] - W.tmpN2[i]
     end
     return dest
 end
-_ldiv!(dest, W, A, B) = _ldiv!(dest, W, lu(A), B)
+_ldiv!(dest, W, A, B) = _ldiv!(dest, W, defaultfactor(W, A), B)
+
+defaultfactor(::AbstractWoodbury, A) = lu(A)
 
 det(W::AbstractWoodbury) = det(W.A)*det(W.C)/det(W.Cp)
 logdet(W::AbstractWoodbury) = logdet(W.A) + logdet(W.C) - logdet(W.Cp)
diff --git a/src/symwoodbury.jl b/src/symwoodbury.jl
index 7fabb78..5dfb142 100644
--- a/src/symwoodbury.jl
+++ b/src/symwoodbury.jl
@@ -91,11 +91,9 @@ inv(O::SymWoodbury{T,AType,BType,DType}) where {T<:Any, AType<:Any, BType<:Any,
 inv(O::SymWoodbury{T,AType,BType,DType}) where {T<:Any, AType<:Any, BType<:Any, DType<:SparseMatrixCSC} =
   calc_inv(O.A, O.B, Matrix(O.D))
 
-_ldiv!(dest, W::SymWoodbury, A::AbstractMatrix, B) = _ldiv!(dest, W, defaultfactor(A), B)
-
-defaultfactor(A) = bunchkaufman(A, true)
-defaultfactor(A::SymTridiagonal) = ldlt(A)
-defaultfactor(A::SparseMatrixCSC) = lu(A)
+defaultfactor(::SymWoodbury, A) = bunchkaufman(A, true)
+defaultfactor(::SymWoodbury, A::SymTridiagonal) = ldlt(A)
+defaultfactor(::SymWoodbury, A::SparseMatrixCSC) = lu(A)
 
 """
     partialInv(O)
diff --git a/test/runtests.jl b/test/runtests.jl
index 0ed9461..9b8af28 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,12 @@
 using WoodburyMatrices
 using LinearAlgebra, SparseArrays, Test
+using Aqua
 using Random: seed!
 
+@testset "Aqua" begin
+    Aqua.test_all(WoodburyMatrices)
+end
+
 # helper function for testing logdet
 function randpsd(sidelength)
     Q = randn(sidelength, sidelength)
diff --git a/test/symwoodbury.jl b/test/symwoodbury.jl
index 09468d6..17e743a 100644
--- a/test/symwoodbury.jl
+++ b/test/symwoodbury.jl
@@ -216,4 +216,15 @@ W = SymWoodbury([randpsd(50) for _ in 1:3]...)
 @test logdet(W) ≈ log(det(W)) ≈ logdet(Array(W))
 @test all(logabsdet(W) .≈ logabsdet(Array(W)))
 
+@testset "Diagonal (issue #35)" begin
+    A = Diagonal( Float64[1,2,3,4] )
+    V = [5  6  7  8;
+         9 10 11 12]
+    U = V'
+    D = Matrix( I, 2, 2 )
+
+    W = SymWoodbury( A, U, D )
+    @test W \ [13, 14, 15, 16] ≈ Matrix(W) \ [13, 14, 15, 16]
+end
+
 end # @testset "SymWoodbury"
diff --git a/test/woodbury.jl b/test/woodbury.jl
index b52218e..b3104bd 100644
--- a/test/woodbury.jl
+++ b/test/woodbury.jl
@@ -65,6 +65,11 @@ for elty in (Float32, Float64, ComplexF32, ComplexF64, Int)
     @test W'*v ≈ F'*v
     @test transpose(W)*v ≈ transpose(F)*v
     @test (W + W) \ v ≈ (2*Matrix(W)) \ v
+    # Test a method used for ambiguity resolution
+    if elty<:Union{Float32, Float64}
+        R = randn(Complex{elty}, n, 2)
+        @test W \ R ≈ Matrix(W) \ R
+    end
 
     # Diagonal matrix for A (lu(A) fails)
     D = Diagonal(d)