From a2062ddb15756cb408743a9e07ea869127f2a031 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreasnoackjensen@gmail.com>
Date: Mon, 3 Nov 2014 09:56:43 -0500
Subject: [PATCH] Fix #8857, i.e. element type promotion in \(SymTridiagonal)

---
 base/linalg/ldlt.jl    | 12 +++++++-----
 base/linalg/lu.jl      |  1 +
 base/linalg/tridiag.jl | 12 +-----------
 3 files changed, 9 insertions(+), 16 deletions(-)

diff --git a/base/linalg/ldlt.jl b/base/linalg/ldlt.jl
index 624ba87363d18..09dd7865cbfe0 100644
--- a/base/linalg/ldlt.jl
+++ b/base/linalg/ldlt.jl
@@ -17,10 +17,9 @@ function ldltfact!{T<:Real}(S::SymTridiagonal{T})
     n = size(S,1)
     d = S.dv
     e = S.ev
-    @inbounds for i = 1:n-1
+    @inbounds @simd for i = 1:n-1
         e[i] /= d[i]
         d[i+1] -= abs2(e[i])*d[i]
-        d[i+1] > 0 || throw(PosDefException(i+1))
     end
     return LDLt{T,SymTridiagonal{T}}(S)
 end
@@ -28,6 +27,9 @@ function ldltfact{T}(M::SymTridiagonal{T})
     S = typeof(zero(T)/one(T))
     return S == T ? ldltfact!(copy(M)) : ldltfact!(convert(SymTridiagonal{S}, M))
 end
+
+factorize(S::SymTridiagonal) = ldltfact(S)
+
 function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
     n, nrhs = size(B, 1), size(B, 2)
     size(S,1) == n || throw(DimensionMismatch(""))
@@ -36,18 +38,18 @@ function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
     @inbounds begin
         for i = 2:n
             li1 = l[i-1]
-            for j = 1:nrhs
+            @simd for j = 1:nrhs
                 B[i,j] -= li1*B[i-1,j]
             end
         end
         dn = d[n]
-        for j = 1:nrhs
+        @simd for j = 1:nrhs
             B[n,j] /= dn
         end
         for i = n-1:-1:1
             di = d[i]
             li = l[i]
-            for j = 1:nrhs
+            @simd for j = 1:nrhs
                 B[i,j] /= di
                 B[i,j] -= li*B[i+1,j]
             end
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index 121e20882cad4..e2b9fbddc20db 100644
--- a/base/linalg/lu.jl
+++ b/base/linalg/lu.jl
@@ -218,6 +218,7 @@ function lufact!{T}(A::Tridiagonal{T}; pivot = true)
     end
     LU{T,Tridiagonal{T}}(A, ipiv, convert(BlasInt, info))
 end
+
 factorize(A::Tridiagonal) = lufact(A)
 
 # See dgtts2.f
diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl
index 1c1d13def358d..3650a984e5fac 100644
--- a/base/linalg/tridiag.jl
+++ b/base/linalg/tridiag.jl
@@ -73,13 +73,7 @@ function A_mul_B!(C::StridedVecOrMat, S::SymTridiagonal, B::StridedVecOrMat)
     return C
 end
 
-## Solver
-function \{T<:BlasFloat}(M::SymTridiagonal{T}, rhs::StridedVecOrMat{T})
-    if stride(rhs, 1) == 1
-        return LAPACK.gtsv!(copy(M.ev), copy(M.dv), copy(M.ev), copy(rhs))
-    end
-    solve(Tridiagonal(M), rhs)  # use the Julia "fallback"
-end
+factorize(S::SymTridiagonal) = ldltfact(S)
 
 #Wrap LAPACK DSTE{GR,BZ} to compute eigenvalues
 eigfact!{T<:BlasFloat}(m::SymTridiagonal{T}) = Eigen(LAPACK.stegr!('V', m.dv, m.ev)...)
@@ -274,7 +268,3 @@ function A_mul_B!(C::AbstractVecOrMat, A::Tridiagonal, B::AbstractVecOrMat)
     end
     C
 end
-
-A_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = A_ldiv_B!(lufact!(A), B)
-At_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = At_ldiv_B!(lufact!(A), B)
-Ac_ldiv_B!(A::Tridiagonal,B::AbstractVecOrMat) = Ac_ldiv_B!(lufact!(A), B)