diff --git a/src/diagonal.jl b/src/diagonal.jl
index 7311f6e1..c919f360 100644
--- a/src/diagonal.jl
+++ b/src/diagonal.jl
@@ -601,10 +601,10 @@ function _rdiv!(B::AbstractVecOrMat, A::AbstractVecOrMat, D::Diagonal)
     if (k = length(dd)) != n
         throw(DimensionMismatch(lazy"left hand side has $n columns but D is $k by $k"))
     end
-    @inbounds for j in 1:n
+    @inbounds for j in axes(A,2)
         ddj = dd[j]
         iszero(ddj) && throw(SingularException(j))
-        for i in 1:m
+        for i in axes(A,1)
             B[i, j] = A[i, j] / ddj
         end
     end
@@ -629,7 +629,17 @@ function ldiv!(B::AbstractVecOrMat, D::Diagonal, A::AbstractVecOrMat)
     (m, n) == (m′, n′) || throw(DimensionMismatch(lazy"expect output to be $m by $n, but got $m′ by $n′"))
     j = findfirst(iszero, D.diag)
     isnothing(j) || throw(SingularException(j))
-    @inbounds for j = 1:n, i = 1:m
+    _ldiv_Diagonal_loop!(B, D, A)
+    B
+end
+function _ldiv_Diagonal_loop!(B::AbstractVecOrMat, D::Diagonal, A::AbstractVecOrMat)
+    dd = D.diag
+    @. B = dd \ A
+    B
+end
+function _ldiv_Diagonal_loop!(B::StridedVecOrMat, D::Diagonal, A::StridedVecOrMat)
+    dd = D.diag
+    @inbounds for j in axes(A,2), i in axes(A,1)
         B[i, j] = dd[i] \ A[i, j]
     end
     B