Merge pull request #8908 from JuliaLang/anj/symtri

Fix #8857, i.e. element type promotion in \(SymTridiagonal)

base/linalg/ldlt.jl

@@ -17,17 +17,19 @@ 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
 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

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

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)