fix order of solve in ldiv

Author: shashi

src/solve.jl

@@ -95,23 +95,24 @@ function simplifying_dot(x,y)
     end
 end
 
-function ldiv(A::LU, b)
-    # unit lower triangular solve first:
+function ldiv(AA, A::LU, b)
     L = A.L
     U = A.U
     m, n = size(L)
     x = Vector{Any}(undef, length(b))
     b = b[A.p]
-    for i=1:n
-        sub = simplifying_dot(b[1:i-1], L[i, 1:i-1])
-        x[i] = _iszero(sub) ? b[i] : b[i] - sub
-    end
 
     for i=n:-1:1
         sub = simplifying_dot(b[i+1:end], U[i,i+1:end])
         den = U[i,i]
+        x[i] = _iszero(sub) ? b[i] : b[i] - sub
+        x[i] = _isone(den) ? x[i] : _isone(-den) ? -x[i] : x[i] / den
+    end
+
+    # unit lower triangular solve first:
+    for i=1:n
+        sub = simplifying_dot(x[1:i-1], L[i, 1:i-1])
         x[i] = _iszero(sub) ? x[i] : x[i] - sub
-        x[i] = _isone(U[i,i]) ? x[i] : _isone(-U[i,i]) ? -x[i] : x[i] / U[i,i]
     end
     x
 end