Merge pull request JuliaLang#11 from JuliaLang/master

update to 0e84b09 [ci skip]

base/linalg/factorization.jl

@@ -250,7 +250,8 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular
 
 inv{T<:BlasFloat}(A::LU{T})=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
 
-cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond{T<:BlasFloat}(A::LU{T}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
 
 ####################
 # QR Factorization #

src/gf.c

@@ -823,6 +823,7 @@ static jl_value_t *lookup_match(jl_value_t *a, jl_value_t *b, jl_tuple_t **penv,
     jl_value_t *ti = jl_type_intersection_matching(a, b, penv, tvars);
     if (ti == (jl_value_t*)jl_bottom_type)
         return ti;
+    JL_GC_PUSH1(&ti);
     assert(jl_is_tuple(*penv));
     jl_value_t **ee = (jl_value_t**)alloca(sizeof(void*) * jl_tuple_len(*penv));
     int n=0;
@@ -854,8 +855,10 @@ static jl_value_t *lookup_match(jl_value_t *a, jl_value_t *b, jl_tuple_t **penv,
                   issue #5254
                 */
                 if (val == (jl_value_t*)jl_bottom_type) {
-                    if (!jl_subtype(a, ti, 0))
+                    if (!jl_subtype(a, ti, 0)) {
+                        JL_GC_POP();
                         return (jl_value_t*)jl_bottom_type;
+                    }
                 }
             }
         }
@@ -865,6 +868,7 @@ static jl_value_t *lookup_match(jl_value_t *a, jl_value_t *b, jl_tuple_t **penv,
         memcpy(en->data, ee, n*sizeof(void*));
         *penv = en;
     }
+    JL_GC_POP();
     return ti;
 }
 

test/linalg1.jl

@@ -20,14 +20,16 @@ areal = randn(n,n)/2
 aimg  = randn(n,n)/2
 breal = randn(n,2)/2
 bimg  = randn(n,2)/2
-for eltya in (Float16, Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
-    for eltyb in (Float16, Float32, Float64, Complex32, Complex64, Complex128, Int)
-        a = eltya == Int ? rand(1:5, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
+for eltya in (Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
+    for eltyb in (Float32, Float64, Complex32, Complex64, Complex128, Int)
+        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
         b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
         asym = a'+a                  # symmetric indefinite
         apd  = a'*a                 # symmetric positive-definite
 
-        ε = max(eps(abs(float(one(eltya)))),eps(abs(float(one(eltyb)))))
+        εa = eps(abs(float(one(eltya))))
+        εb = eps(abs(float(one(eltyb))))
+        ε = max(εa,εb)
 
 debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
 
@@ -50,8 +52,8 @@ debug && println("(Automatic) upper Cholesky factor")
         @test norm(apd*x-b)/norm(b) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
 
         @test_approx_eq apd * inv(capd) eye(n)
-        @test_approx_eq_eps a*(capd\(a'*b)) b 8000ε
-        @test_approx_eq det(capd) det(apd)
+        @test norm(a*(capd\(a'*b)) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
+        @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit
         @test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow
 
 debug && println("lower Cholesky factor")
@@ -64,7 +66,7 @@ debug && println("pivoted Choleksy decomposition")
         cpapd = cholfact(apd, pivot=true)
         @test rank(cpapd) == n
         @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
-        @test_approx_eq_eps b apd * (cpapd\b) 15000ε
+        @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
         if isreal(apd)
             @test_approx_eq apd * inv(cpapd) eye(n)
         end
@@ -83,12 +85,13 @@ debug && println("Bunch-Kaufman factors of a pos-def matrix")
     end
 
 debug && println("(Automatic) Square LU decomposition")
+    κ     = cond(a,1)
     lua   = factorize(a)
     l,u,p = lua[:L], lua[:U], lua[:p]
     @test_approx_eq l*u a[p,:]
     @test_approx_eq l[invperm(p),:]*u a
     @test_approx_eq a * inv(lua) eye(n)
-    @test_approx_eq_eps a*(lua\b) b 80ε
+    @test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
 
 debug && println("Thin LU")
     lua   = lufact(a[:,1:5])
@@ -188,10 +191,11 @@ debug && println("Generalized svd")
     end
 
 debug && println("Solve square general system of equations")
+    κ = cond(a,1)
     x = a \ b
-    @test_approx_eq_eps a*x b 80ε
     @test_throws b'\b
     @test_throws b\b'
+    @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit!
 
 debug && println("Solve upper triangular system")
     x = triu(a) \ b

test/linalg2.jl

@@ -43,6 +43,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
         V = convert(Matrix{elty}, V)
         C = convert(Matrix{elty}, C)
     end
+    ε = eps(abs2(float(one(elty))))
     T = Tridiagonal(dl, d, du)
     @test size(T, 1) == n
     @test size(T) == (n, n)
@@ -95,8 +96,8 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
     F = full(W)
     @test_approx_eq W*v F*v
     iFv = F\v
-    @test_approx_eq W\v iFv
-    @test_approx_eq det(W) det(F)
+    @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong
+    @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
     iWv = similar(iFv)
     if elty != Int
         Base.LinAlg.solve!(iWv, W, v)
@@ -208,7 +209,7 @@ end
 
 #Triangular matrices
 n=12
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     A = convert(Matrix{elty}, randn(n, n))
     b = convert(Matrix{elty}, randn(n, 2))
     if elty <: Complex
@@ -245,7 +246,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt
 end
 
 #Tridiagonal matrices
-for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty})
+for relty in (Float32, Float64), elty in (relty, Complex{relty})
     a = convert(Vector{elty}, randn(n-1))
     b = convert(Vector{elty}, randn(n))
     c = convert(Vector{elty}, randn(n-1))
@@ -263,11 +264,10 @@ for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty})
 end
 
 #SymTridiagonal (symmetric tridiagonal) matrices
-for relty in (Float16, Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work
+for relty in (Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work
     a = convert(Vector{elty}, randn(n))
     b = convert(Vector{elty}, randn(n-1))
     if elty <: Complex
-        relty==Float16 && continue
         a += im*convert(Vector{elty}, randn(n))
         b += im*convert(Vector{elty}, randn(n-1))
     end
@@ -306,7 +306,7 @@ for elty in (Float32, Float64)
 end
 
 #Bidiagonal matrices
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     dv = convert(Vector{elty}, randn(n))
     ev = convert(Vector{elty}, randn(n-1))
     b = convert(Matrix{elty}, randn(n, 2))
@@ -361,7 +361,7 @@ end
 
 #Diagonal matrices
 n=12
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
     d=convert(Vector{elty}, randn(n))
     v=convert(Vector{elty}, randn(n))
     U=convert(Matrix{elty}, randn(n,n))
@@ -529,7 +529,7 @@ end
 nnorm = 1000
 mmat = 100
 nmat = 80
-for elty in (Float16, Float32, Float64, BigFloat, Complex{Float16}, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
+for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
     debug && println(elty)
 
     ## Vector