Merge pull request #10215 from JuliaLang/anj/lu

Some changes to LU and triangular solves

base/linalg/lu.jl

@@ -16,7 +16,7 @@ function lufact!{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Union(Type{Val{false}
     return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
 end
 lufact!(A::StridedMatrix, pivot::Union(Type{Val{false}}, Type{Val{true}}) = Val{true}) = generic_lufact!(A, pivot)
-function generic_lufact!{T}(A::StridedMatrix{T}, pivot::Union(Type{Val{false}}, Type{Val{true}}) = Val{true})
+function generic_lufact!{T,Pivot}(A::StridedMatrix{T}, ::Type{Val{Pivot}} = Val{true})
     m, n = size(A)
     minmn = min(m,n)
     info = 0
@@ -25,7 +25,7 @@ function generic_lufact!{T}(A::StridedMatrix{T}, pivot::Union(Type{Val{false}},
         for k = 1:minmn
             # find index max
             kp = k
-            if pivot==Val{true}
+            if Pivot
                 amax = real(zero(T))
                 for i = k:m
                     absi = abs(A[i,k])
@@ -63,8 +63,26 @@ function generic_lufact!{T}(A::StridedMatrix{T}, pivot::Union(Type{Val{false}},
     end
     LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
 end
-lufact{T<:BlasFloat}(A::AbstractMatrix{T}, pivot::Union(Type{Val{false}}, Type{Val{true}}) = Val{true}) = lufact!(copy(A), pivot)
-lufact{T}(A::AbstractMatrix{T}, pivot::Union(Type{Val{false}}, Type{Val{true}}) = Val{true}) = (S = typeof(zero(T)/one(T)); S != T ? lufact!(convert(AbstractMatrix{S}, A), pivot) : lufact!(copy(A), pivot))
+
+# floating point types doesn't have to be promoted for LU, but should default to pivoting
+lufact{T<:FloatingPoint}(A::Union(AbstractMatrix{T},AbstractMatrix{Complex{T}}), pivot::Union(Type{Val{false}}, Type{Val{true}}) = Val{true}) = lufact!(copy(A), pivot)
+
+# for all other types we must promote to a type which is stable under division
+function lufact{T}(A::AbstractMatrix{T}, pivot::Union(Type{Val{false}}, Type{Val{true}}))
+    S = typeof(zero(T)/one(T))
+    lufact!(copy_oftype(A, S), pivot)
+end
+# We can't assume an ordered field so we first try without pivoting
+function lufact{T}(A::AbstractMatrix{T})
+    S = typeof(zero(T)/one(T))
+    F = lufact!(copy_oftype(A, S), Val{false})
+    if F.info == 0
+        return F
+    else
+        return lufact!(copy_oftype(A, S), Val{true})
+    end
+end
+
 lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt))
 lufact(F::LU) = F
 

base/linalg/triangular.jl

@@ -560,48 +560,60 @@ end
 
 #Generic solver using naive substitution
 function naivesub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector=b)
-    N = size(A, 2)
-    N == length(b) == length(x) || throw(DimensionMismatch())
-    for j = N:-1:1
-        x[j] = b[j]
-        for k = j+1:1:N
-            x[j] -= A[j,k] * x[k]
+    n = size(A, 2)
+    n == length(b) == length(x) || throw(DimensionMismatch())
+    for j = n:-1:1
+        xj = b[j]
+        for k = j+1:1:n
+            xj -= A[j,k] * x[k]
+        end
+        Ajj = A[j,j]
+        if Ajj == zero(Ajj)
+            throw(SingularException(j))
+        else
+            x[j] = Ajj\xj
         end
-        x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
     end
     x
 end
 function naivesub!(A::UnitUpperTriangular, b::AbstractVector, x::AbstractVector=b)
-    N = size(A, 2)
-    N == length(b) == length(x) || throw(DimensionMismatch())
-    for j = N:-1:1
-        x[j] = b[j]
-        for k = j+1:1:N
-            x[j] -= A[j,k] * x[k]
+    n = size(A, 2)
+    n == length(b) == length(x) || throw(DimensionMismatch())
+    for j = n:-1:1
+        xj = b[j]
+        for k = j+1:1:n
+            xj -= A[j,k] * x[k]
         end
+        x[j] = xj
     end
     x
 end
 function naivesub!(A::LowerTriangular, b::AbstractVector, x::AbstractVector=b)
-    N = size(A, 2)
-    N == length(b) == length(x) || throw(DimensionMismatch())
-    for j = 1:N
-        x[j] = b[j]
+    n = size(A, 2)
+    n == length(b) == length(x) || throw(DimensionMismatch())
+    for j = 1:n
+        xj = b[j]
         for k = 1:j-1
-            x[j] -= A[j,k] * x[k]
+            xj -= A[j,k] * x[k]
+        end
+        Ajj = A[j,j]
+        if Ajj == zero(Ajj)
+            throw(SingularException(j))
+        else
+            x[j] = Ajj\xj
         end
-        x[j] = A[j,j]==0 ? throw(SingularException(j)) : A[j,j]\x[j]
     end
     x
 end
 function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector=b)
-    N = size(A, 2)
-    N == length(b) == length(x) || throw(DimensionMismatch())
-    for j = 1:N
-        x[j] = b[j]
+    n = size(A, 2)
+    n == length(b) == length(x) || throw(DimensionMismatch())
+    for j = 1:n
+        xj = b[j]
         for k = 1:j-1
-            x[j] -= A[j,k] * x[k]
+            xj -= A[j,k] * x[k]
         end
+        x[j] = xj
     end
     x
 end