From 0eaa06d28a4193a3fd935877f939ea64eab11c8a Mon Sep 17 00:00:00 2001 From: kshyatt Date: Fri, 31 Jul 2015 17:51:59 -0700 Subject: [PATCH] Moved linalg tests into appropriately named files --- test/blas.jl | 33 +++ test/choosetests.jl | 13 +- test/linalg/bunchkaufman.jl | 55 +++++ test/linalg/cholesky.jl | 23 +- test/linalg/dense.jl | 434 ++++++++++++++++++++++++++++++++++++ test/linalg/eigen.jl | 62 ++++++ test/linalg/generic.jl | 77 +++++++ test/linalg/lapack.jl | 25 +++ test/linalg/matmul.jl | 134 +++++++++++ test/linalg/qr.jl | 126 +++++++++++ test/linalg/schur.jl | 70 ++++++ test/linalg/special.jl | 105 +++++++++ test/linalg/svd.jl | 50 +++++ test/linalg/tridiag.jl | 114 ++++++++++ test/linalg1.jl | 323 --------------------------- test/linalg2.jl | 395 -------------------------------- test/linalg3.jl | 333 --------------------------- test/linalg4.jl | 149 ------------- test/sparsedir/sparse.jl | 16 ++ 19 files changed, 1321 insertions(+), 1216 deletions(-) create mode 100644 test/linalg/bunchkaufman.jl create mode 100644 test/linalg/eigen.jl create mode 100644 test/linalg/matmul.jl create mode 100644 test/linalg/qr.jl create mode 100644 test/linalg/schur.jl create mode 100644 test/linalg/special.jl create mode 100644 test/linalg/svd.jl delete mode 100644 test/linalg1.jl delete mode 100644 test/linalg2.jl delete mode 100644 test/linalg3.jl delete mode 100644 test/linalg4.jl diff --git a/test/blas.jl b/test/blas.jl index 399d8fb90d973..a029d74b2af36 100644 --- a/test/blas.jl +++ b/test/blas.jl @@ -1,6 +1,39 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license import Base.LinAlg + +srand(100) +# syr2k! and her2k! +for elty in (Float32, Float64, Complex64, Complex128) + U = randn(5,2) + V = randn(5,2) + if elty == Complex64 || elty == Complex128 + U = complex(U, U) + V = complex(V, V) + end + U = convert(Array{elty, 2}, U) + V = convert(Array{elty, 2}, V) + @test_approx_eq tril(LinAlg.BLAS.syr2k('L','N',U,V)) tril(U*V.' + V*U.') + @test_approx_eq triu(LinAlg.BLAS.syr2k('U','N',U,V)) triu(U*V.' + V*U.') + @test_approx_eq tril(LinAlg.BLAS.syr2k('L','T',U,V)) tril(U.'*V + V.'*U) + @test_approx_eq triu(LinAlg.BLAS.syr2k('U','T',U,V)) triu(U.'*V + V.'*U) +end + +for elty in (Complex64, Complex128) + U = randn(5,2) + V = randn(5,2) + if elty == Complex64 || elty == Complex128 + U = complex(U, U) + V = complex(V, V) + end + U = convert(Array{elty, 2}, U) + V = convert(Array{elty, 2}, V) + @test_approx_eq tril(LinAlg.BLAS.her2k('L','N',U,V)) tril(U*V' + V*U') + @test_approx_eq triu(LinAlg.BLAS.her2k('U','N',U,V)) triu(U*V' + V*U') + @test_approx_eq tril(LinAlg.BLAS.her2k('L','C',U,V)) tril(U'*V + V'*U) + @test_approx_eq triu(LinAlg.BLAS.her2k('U','C',U,V)) triu(U'*V + V'*U) +end + ## BLAS tests - testing the interface code to BLAS routines for elty in [Float32, Float64, Complex64, Complex128] diff --git a/test/choosetests.jl b/test/choosetests.jl index 7bf4a49aa82a9..f85a553e79d4a 100644 --- a/test/choosetests.jl +++ b/test/choosetests.jl @@ -62,12 +62,13 @@ function choosetests(choices = []) tests = testnames end - linalgtests = ["linalg1", "linalg2", "linalg3", "linalg4", - "linalg/lapack", "linalg/triangular", "linalg/tridiag", - "linalg/bidiag", "linalg/diagonal", "linalg/dense", - "linalg/pinv", "linalg/givens", "linalg/cholesky", - "linalg/lu", "linalg/symmetric", "linalg/generic", - "linalg/uniformscaling"] + linalgtests = ["linalg/triangular", "linalg/qr", "linalg/dense", + "linalg/matmul", "linalg/schur", "linalg/special", + "linalg/eigen", "linalg/bunchkaufman", "linalg/svd", + "linalg/lapack", "linalg/tridiag", "linalg/bidiag", + "linalg/diagonal", "linalg/pinv", "linalg/givens", + "linalg/cholesky", "linalg/lu", "linalg/symmetric", + "linalg/generic", "linalg/uniformscaling"] if Base.USE_GPL_LIBS push!(linalgtests, "linalg/arnoldi") end diff --git a/test/linalg/bunchkaufman.jl b/test/linalg/bunchkaufman.jl new file mode 100644 index 0000000000000..64a5db6899b50 --- /dev/null +++ b/test/linalg/bunchkaufman.jl @@ -0,0 +1,55 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +debug = false +using Base.Test + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 +a2real = randn(n,n)/2 +a2img = randn(n,n)/2 +breal = randn(n,2)/2 +bimg = randn(n,2)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + + for eltyb in (Float32, Float64, Complex64, Complex128, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + +debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") + +debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix") + bc1 = factorize(asym) + @test_approx_eq inv(bc1) * asym eye(n) + @test_approx_eq_eps asym * (bc1\b) b 1000ε + @test_approx_eq inv(a.' + a) * (a.' + a) eye(n) + @test size(bc1) == size(bc1.LD) + @test size(bc1,1) == size(bc1.LD,1) + @test size(bc1,2) == size(bc1.LD,2) + if eltya <: BlasReal + @test_throws ArgumentError bkfact(a) + end +debug && println("Bunch-Kaufman factors of a pos-def matrix") + bc2 = bkfact(apd) + @test_approx_eq inv(bc2) * apd eye(n) + @test_approx_eq_eps apd * (bc2\b) b 150000ε + @test ishermitian(bc2) == !issym(bc2) + + end +end diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl index 72013e4ecaff2..57559292997f7 100644 --- a/test/linalg/cholesky.jl +++ b/test/linalg/cholesky.jl @@ -14,16 +14,6 @@ n2 = 2*n1 srand(1234321) -a = rand(n,n) -for elty in (Float32, Float64, Complex64, Complex128) - a = convert(Matrix{elty}, a) - # cond - @test_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01 - @test_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01 - @test_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01 - @test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01 -end - areal = randn(n,n)/2 aimg = randn(n,n)/2 a2real = randn(n,n)/2 @@ -112,3 +102,16 @@ begin @test sum(sum(norm, L*L' - XX)) < eps() @test sum(sum(norm, U'*U - XX)) < eps() end + +# Test generic cholfact! +for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) + if elty <: Complex + A = complex(randn(5,5), randn(5,5)) + else + A = randn(5,5) + end + A = convert(Matrix{elty}, A'A) + for ul in (:U, :L) + @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{Val{ul}}},copy(A), Val{ul})) + end +end diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index ca9aed280c307..0c062ed792291 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -1,6 +1,440 @@ # This file is a part of Julia. License is MIT: http://julialang.org/license +debug = false using Base.Test # Check that non-floats are correctly promoted @test_approx_eq [1 0 0; 0 1 0]\[1,1] [1;1;0] + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +a = rand(n,n) +for elty in (Float32, Float64, Complex64, Complex128) + a = convert(Matrix{elty}, a) + # cond + @test_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01 + @test_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01 + @test_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01 + @test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01 +end + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 +a2real = randn(n,n)/2 +a2img = randn(n,n)/2 +breal = randn(n,2)/2 +bimg = randn(n,2)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + + for eltyb in (Float32, Float64, Complex64, Complex128, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + +debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") + +debug && println("Solve square general system of equations") + κ = cond(a,1) + x = a \ b + @test_throws DimensionMismatch b'\b + @test_throws DimensionMismatch b\b' + @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit! + + +debug && println("Test nullspace") + a15null = nullspace(a[:,1:n1]') + @test rank([a[:,1:n1] a15null]) == 10 + @test_approx_eq_eps norm(a[:,1:n1]'a15null, Inf) zero(eltya) 300ε + @test_approx_eq_eps norm(a15null'a[:,1:n1], Inf) zero(eltya) 400ε + @test size(nullspace(b), 2) == 0 + end # for eltyb + +debug && println("Test pinv") + pinva15 = pinv(a[:,1:n1]) + @test_approx_eq a[:,1:n1]*pinva15*a[:,1:n1] a[:,1:n1] + @test_approx_eq pinva15*a[:,1:n1]*pinva15 pinva15 + +debug && println("Lyapunov/Sylvester") + let + x = lyap(a, a2) + @test_approx_eq -a2 a*x + x*a' + x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n]) + @test_approx_eq -a2[1:3, 4:n] a[1:3, 1:3]*x2 + x2*a[4:n, 4:n] + end + +debug && println("Matrix square root") + asq = sqrtm(a) + @test_approx_eq asq*asq a + asymsq = sqrtm(asym) + @test_approx_eq asymsq*asymsq asym +end # for eltya + +# test triu/tril bounds checking +A = rand(5,7) +@test_throws(BoundsError,triu(A,8)) +@test_throws(BoundsError,triu(A,-6)) +@test_throws(BoundsError,tril(A,8)) +@test_throws(BoundsError,tril(A,-6)) + +# Test gradient +for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128) + if elty <: Real + x = convert(Vector{elty}, [1:3;]) + g = ones(elty, 3) + else + x = convert(Vector{elty}, complex([1:3;], [1:3;])) + g = convert(Vector{elty}, complex(ones(3), ones(3))) + end + @test_approx_eq gradient(x) g +end + +# Tests norms +nnorm = 10 +mmat = 10 +nmat = 8 +for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt) + debug && println(elty) + + ## Vector + x = ones(elty,10) + xs = sub(x,1:2:10) + @test_approx_eq norm(x, -Inf) 1 + @test_approx_eq norm(x, -1) 1/10 + @test_approx_eq norm(x, 0) 10 + @test_approx_eq norm(x, 1) 10 + @test_approx_eq norm(x, 2) sqrt(10) + @test_approx_eq norm(x, 3) cbrt(10) + @test_approx_eq norm(x, Inf) 1 + @test_approx_eq norm(xs, -Inf) 1 + @test_approx_eq norm(xs, -1) 1/5 + @test_approx_eq norm(xs, 0) 5 + @test_approx_eq norm(xs, 1) 5 + @test_approx_eq norm(xs, 2) sqrt(5) + @test_approx_eq norm(xs, 3) cbrt(5) + @test_approx_eq norm(xs, Inf) 1 + + ## Number + norm(x[1:1]) === norm(x[1], -Inf) + norm(x[1:1]) === norm(x[1], 0) + norm(x[1:1]) === norm(x[1], 1) + norm(x[1:1]) === norm(x[1], 2) + norm(x[1:1]) === norm(x[1], Inf) + + for i = 1:10 + x = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : + elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : + convert(Vector{elty}, randn(nnorm)) + xs = sub(x,1:2:nnorm) + y = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : + elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : + convert(Vector{elty}, randn(nnorm)) + ys = sub(y,1:2:nnorm) + α = elty <: Integer ? randn() : + elty <: Complex ? convert(elty, complex(randn(),randn())) : + convert(elty, randn()) + # Absolute homogeneity + @test_approx_eq norm(α*x,-Inf) abs(α)*norm(x,-Inf) + @test_approx_eq norm(α*x,-1) abs(α)*norm(x,-1) + @test_approx_eq norm(α*x,1) abs(α)*norm(x,1) + @test_approx_eq norm(α*x) abs(α)*norm(x) # two is default + @test_approx_eq norm(α*x,3) abs(α)*norm(x,3) + @test_approx_eq norm(α*x,Inf) abs(α)*norm(x,Inf) + + @test_approx_eq norm(α*xs,-Inf) abs(α)*norm(xs,-Inf) + @test_approx_eq norm(α*xs,-1) abs(α)*norm(xs,-1) + @test_approx_eq norm(α*xs,1) abs(α)*norm(xs,1) + @test_approx_eq norm(α*xs) abs(α)*norm(xs) # two is default + @test_approx_eq norm(α*xs,3) abs(α)*norm(xs,3) + @test_approx_eq norm(α*xs,Inf) abs(α)*norm(xs,Inf) + + # Triangle inequality + @test norm(x + y,1) <= norm(x,1) + norm(y,1) + @test norm(x + y) <= norm(x) + norm(y) # two is default + @test norm(x + y,3) <= norm(x,3) + norm(y,3) + @test norm(x + y,Inf) <= norm(x,Inf) + norm(y,Inf) + + @test norm(xs + ys,1) <= norm(xs,1) + norm(ys,1) + @test norm(xs + ys) <= norm(xs) + norm(ys) # two is default + @test norm(xs + ys,3) <= norm(xs,3) + norm(ys,3) + @test norm(xs + ys,Inf) <= norm(xs,Inf) + norm(ys,Inf) + + # Against vectorized versions + @test_approx_eq norm(x,-Inf) minimum(abs(x)) + @test_approx_eq norm(x,-1) inv(sum(1./abs(x))) + @test_approx_eq norm(x,0) sum(x .!= 0) + @test_approx_eq norm(x,1) sum(abs(x)) + @test_approx_eq norm(x) sqrt(sum(abs2(x))) + @test_approx_eq norm(x,3) cbrt(sum(abs(x).^3.)) + @test_approx_eq norm(x,Inf) maximum(abs(x)) + end + ## Matrix (Operator) + A = ones(elty,10,10) + As = sub(A,1:5,1:5) + @test_approx_eq norm(A, 1) 10 + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(A, 2) 10 + @test_approx_eq norm(A, Inf) 10 + @test_approx_eq norm(As, 1) 5 + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(As, 2) 5 + @test_approx_eq norm(As, Inf) 5 + + for i = 1:10 + A = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : + elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : + convert(Matrix{elty}, randn(mmat, nmat)) + As = sub(A,1:nmat,1:nmat) + B = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : + elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : + convert(Matrix{elty}, randn(mmat, nmat)) + Bs = sub(B,1:nmat,1:nmat) + α = elty <: Integer ? randn() : + elty <: Complex ? convert(elty, complex(randn(),randn())) : + convert(elty, randn()) + + # Absolute homogeneity + @test_approx_eq norm(α*A,1) abs(α)*norm(A,1) + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*A) abs(α)*norm(A) # two is default + @test_approx_eq norm(α*A,Inf) abs(α)*norm(A,Inf) + + @test_approx_eq norm(α*As,1) abs(α)*norm(As,1) + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*As) abs(α)*norm(As) # two is default + @test_approx_eq norm(α*As,Inf) abs(α)*norm(As,Inf) + + # Triangle inequality + @test norm(A + B,1) <= norm(A,1) + norm(B,1) + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A + B) <= norm(A) + norm(B) # two is default + @test norm(A + B,Inf) <= norm(A,Inf) + norm(B,Inf) + + @test norm(As + Bs,1) <= norm(As,1) + norm(Bs,1) + elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(As + Bs) <= norm(As) + norm(Bs) # two is default + @test norm(As + Bs,Inf) <= norm(As,Inf) + norm(Bs,Inf) + + # vecnorm: + for p = -2:3 + @test norm(reshape(A, length(A)), p) == vecnorm(A, p) + end + + # issue #10234 + if elty <: AbstractFloat || elty <: Complex + let z = zeros(elty, 100) + z[1] = -Inf + for p in [-2,-1.5,-1,-0.5,0.5,1,1.5,2,Inf] + @test norm(z, p) == (p < 0 ? 0 : Inf) + @test norm(elty[Inf],p) == Inf + end + end + end + end +end + +# issue #10234 +@test norm(Any[Inf],-2) == norm(Any[Inf],-1) == norm(Any[Inf],1) == norm(Any[Inf],1.5) == norm(Any[Inf],2) == norm(Any[Inf],Inf) == Inf + +# overflow/underflow in norms: +@test_approx_eq norm(Float64[1e-300, 1], -3)*1e300 1 +@test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1 + +## Issue related tests +# issue #1447 +let + A = [1.+0.im 0; 0 1] + B = pinv(A) + for i = 1:4 + @test_approx_eq A[i] B[i] + end +end + +# issue #2246 +let + A = [1 2 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0] + Asq = sqrtm(A) + @test_approx_eq Asq*Asq A + A2 = sub(A, 1:2, 1:2) + A2sq = sqrtm(A2) + @test_approx_eq A2sq*A2sq A2 +end + +let + N = 3 + @test_approx_eq log(det(eye(N))) logdet(eye(N)) +end + +# issue #2637 +let + a = [1, 2, 3] + b = [4, 5, 6] + @test kron(eye(2),eye(2)) == eye(4) + @test kron(a,b) == [4,5,6,8,10,12,12,15,18] + @test kron(a',b') == [4 5 6 8 10 12 12 15 18] + @test kron(a,b') == [4 5 6; 8 10 12; 12 15 18] + @test kron(a',b) == [4 8 12; 5 10 15; 6 12 18] + @test kron(a,eye(2)) == [1 0; 0 1; 2 0; 0 2; 3 0; 0 3] + @test kron(eye(2),a) == [ 1 0; 2 0; 3 0; 0 1; 0 2; 0 3] + @test kron(eye(2),2) == 2*eye(2) + @test kron(3,eye(3)) == 3*eye(3) + @test kron(a,2) == [2, 4, 6] + @test kron(b',2) == [8 10 12] +end + +# issue #4796 +let + dim=2 + S=zeros(Complex,dim,dim) + T=zeros(Complex,dim,dim) + T[:] = 1 + z = 2.5 + 1.5im + S[1] = z + @test S*T == [z z; 0 0] +end + +# Matrix exponential +for elty in (Float32, Float64, Complex64, Complex128) + A1 = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4]) + eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181 127.781085523182; + 183.765138646367 183.765138646366 163.679601723179; + 71.797032399996 91.8825693231832 111.968106246371]') + @test_approx_eq expm(A1) eA1 + + A2 = convert(Matrix{elty}, + [29.87942128909879 0.7815750847907159 -2.289519314033932; + 0.7815750847907159 25.72656945571064 8.680737820540137; + -2.289519314033932 8.680737820540137 34.39400925519054]) + eA2 = convert(Matrix{elty}, + [ 5496313853692458.0 -18231880972009236.0 -30475770808580460.0; + -18231880972009252.0 60605228702221920.0 101291842930249760.0; + -30475770808580480.0 101291842930249728.0 169294411240851968.0]) + @test_approx_eq expm(A2) eA2 + + A3 = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52]) + eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812 -4.934938326092; + 0.367879439109187 1.47151775849686 1.10363831732856; + 0.135335281175235 0.406005843524598 0.541341126763207]') + @test_approx_eq expm(A3) eA3 + + # Hessenberg + @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty}, + [4.000000000000000 -1.414213562373094 -1.414213562373095 + -1.414213562373095 4.999999999999996 -0.000000000000000 + 0 -0.000000000000002 3.000000000000000]) + @test_throws KeyError hessfact(A1)[:Z] + @test_approx_eq full(hessfact(A1)) A1 + @test_approx_eq full(Base.LinAlg.HessenbergQ(hessfact(A1))) full(hessfact(A1)[:Q]) +end + +for elty in (Float64, Complex{Float64}) + + A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps(); + 1/3 1/4 1/5 1/6; + 1/4 1/5 1/6 1/7; + 1/5 1/6 1/7 1/8]) + @test_approx_eq expm(logm(A4)) A4 + + OLD_STDERR = STDERR + rd,wr = redirect_stderr() + A5 = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1]) + @test_approx_eq expm(logm(A5)) A5 + s = readline(rd) + @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.") + + A6 = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11]) + @test_approx_eq expm(logm(A6)) A6 + s = readline(rd) + @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.") + redirect_stderr(OLD_STDERR) + + A7 = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]) + @test_approx_eq expm(logm(A7)) A7 + +end + +A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]; +@test_approx_eq expm(logm(A8)) A8 + +# issue 5116 +A9 = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0] +eA9 = [-0.999786072879326 -0.065407069689389 0.0 0.0 + 0.006540706968939 -0.999786072879326 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.013081413937878 -3.999572145758650 0.0 1.0] +@test_approx_eq expm(A9) eA9 + +# issue 5116 +A10 = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.] +eA10 = [ 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im + 0.0+0.0im 1.543080634815244+0.0im 0.0-1.175201193643801im 0.0+0.0im + 0.0+0.0im 0.0+1.175201193643801im 1.543080634815243+0.0im 0.0+0.0im + 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im] +@test_approx_eq expm(A10) eA10 + +# issue #7181 +A = [ 1 5 9 + 2 6 10 + 3 7 11 + 4 8 12 ] +@test_throws BoundsError diag(A, -5) +@test diag(A,-4) == [] +@test diag(A,-3) == [4] +@test diag(A,-2) == [3,8] +@test diag(A,-1) == [2,7,12] +@test diag(A, 0) == [1,6,11] +@test diag(A, 1) == [5,10] +@test diag(A, 2) == [9] +@test diag(A, 3) == [] +@test_throws BoundsError diag(A, 4) + +@test diag(zeros(0,0)) == [] +@test_throws BoundsError diag(zeros(0,0),1) +@test_throws BoundsError diag(zeros(0,0),-1) + +@test diag(zeros(1,0)) == [] +@test diag(zeros(1,0),-1) == [] +@test_throws BoundsError diag(zeros(1,0),1) +@test_throws BoundsError diag(zeros(1,0),-2) + +@test diag(zeros(0,1)) == [] +@test diag(zeros(0,1),1) == [] +@test_throws BoundsError diag(zeros(0,1),-1) +@test_throws BoundsError diag(zeros(0,1),2) + +## Least squares solutions +a = [ones(20) 1:20 1:20] +b = reshape(eye(8, 5), 20, 2) +for elty in (Float32, Float64, Complex64, Complex128) + a = convert(Matrix{elty}, a) + b = convert(Matrix{elty}, b) + + # Vector rhs + x = a[:,1:2]\b[:,1] + @test_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384) + + # Matrix rhs + x = a[:,1:2]\b + @test_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031) + + # Rank deficient + x = a\b + @test_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031) + + # Underdetermined minimum norm + x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) + @test_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5]) + + # symmetric, positive definite + @test_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3]) + + # symmetric, indefinite + @test_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3) +end diff --git a/test/linalg/eigen.jl b/test/linalg/eigen.jl new file mode 100644 index 0000000000000..e2156399de62f --- /dev/null +++ b/test/linalg/eigen.jl @@ -0,0 +1,62 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +debug = false +using Base.Test + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + +debug && println("\ntype of a: ", eltya,"\n") +debug && println("non-symmetric eigen decomposition") + d,v = eig(a) + for i in 1:size(a,2) @test_approx_eq a*v[:,i] d[i]*v[:,i] end + f = eigfact(a) + @test_approx_eq det(a) det(f) + @test_approx_eq inv(a) inv(f) + + num_fact = eigfact(one(eltya)) + @test num_fact.values[1] == one(eltya) + h = a + a' + @test_approx_eq minimum(eigvals(h)) eigmin(h) + @test_approx_eq maximum(eigvals(h)) eigmax(h) + @test_throws DomainError eigmin(a - a') + @test_throws DomainError eigmax(a - a') + +debug && println("symmetric generalized eigenproblem") + asym_sg = asym[1:n1, 1:n1] + a_sg = a[:,n1+1:n2] + f = eigfact(asym_sg, a_sg'a_sg) + eig(asym_sg, a_sg'a_sg) # same result, but checks that method works + @test_approx_eq asym_sg*f[:vectors] scale(a_sg'a_sg*f[:vectors], f[:values]) + @test_approx_eq f[:values] eigvals(asym_sg, a_sg'a_sg) + @test_approx_eq_eps prod(f[:values]) prod(eigvals(asym_sg/(a_sg'a_sg))) 200ε + @test eigvecs(asym_sg, a_sg'a_sg) == f[:vectors] + @test_throws KeyError f[:Z] + +debug && println("Non-symmetric generalized eigenproblem") + a1_nsg = a[1:n1, 1:n1] + a2_nsg = a[n1+1:n2, n1+1:n2] + f = eigfact(a1_nsg, a2_nsg) + eig(a1_nsg, a2_nsg) # same result, but checks that method works + @test_approx_eq a1_nsg*f[:vectors] scale(a2_nsg*f[:vectors], f[:values]) + @test_approx_eq f[:values] eigvals(a1_nsg, a2_nsg) + @test_approx_eq_eps prod(f[:values]) prod(eigvals(a1_nsg/a2_nsg)) 50000ε + @test eigvecs(a1_nsg, a2_nsg) == f[:vectors] + @test_throws KeyError f[:Z] +end diff --git a/test/linalg/generic.jl b/test/linalg/generic.jl index a33a90013e57a..329f17caeb967 100644 --- a/test/linalg/generic.jl +++ b/test/linalg/generic.jl @@ -29,3 +29,80 @@ for elty in (Int, Rational{BigInt}, Float32, Float64, BigFloat, Complex{Float32} @test logabsdet(convert(Matrix{elty}, -eye(n)))[2] == -1 end end + +# test diff, throw ArgumentError for invalid dimension argument +let X = [3 9 5; + 7 4 2; + 2 1 10] + @test diff(X,1) == [4 -5 -3; -5 -3 8] + @test diff(X,2) == [6 -4; -3 -2; -1 9] + @test_throws ArgumentError diff(X,3) + @test_throws ArgumentError diff(X,-1) +end + +x = float([1:12;]) +y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99] +@test_approx_eq linreg(x,y) [2.5559090909090867, 1.6960139860139862] + +# test diag +let A = eye(4) + @test diag(A) == ones(4) + @test diag(sub(A, 1:3, 1:3)) == ones(3) +end + +# test generic axpy +x = ['a','b','c','d','e'] +y = ['a','b','c','d','e'] +α = 'f' +@test_throws DimensionMismatch Base.LinAlg.axpy!(α,x,['g']) +@test_throws BoundsError Base.LinAlg.axpy!(α,x,collect(-1:5),y,collect(1:7)) +@test_throws BoundsError Base.LinAlg.axpy!(α,x,collect(1:7),y,collect(1:7)) + + + + +# 2-argument version of scale +a = reshape([1.:6;], (2,3)) +@test scale(a, 5.) == a*5 +@test scale(5., a) == a*5 +@test scale(a, [1.; 2.; 3.]) == a.*[1 2 3] +@test scale([1.; 2.], a) == a.*[1; 2] +@test scale(a, [1; 2; 3]) == a.*[1 2 3] +@test scale([1; 2], a) == a.*[1; 2] +@test_throws DimensionMismatch scale(a, ones(2)) +@test_throws DimensionMismatch scale(ones(3), a) + +# 2-argument version of scale! +@test scale!(copy(a), 5.) == a*5 +@test scale!(5., copy(a)) == a*5 +b = randn(Base.LinAlg.SCAL_CUTOFF) # make sure we try BLAS path +@test scale!(copy(b), 5.) == b*5 +@test scale!(copy(a), [1.; 2.; 3.]) == a.*[1 2 3] +@test scale!([1.; 2.], copy(a)) == a.*[1; 2] +@test scale!(copy(a), [1; 2; 3]) == a.*[1 2 3] +@test scale!([1; 2], copy(a)) == a.*[1; 2] +@test_throws DimensionMismatch scale!(a, ones(2)) +@test_throws DimensionMismatch scale!(ones(3), a) + +# 3-argument version of scale! +@test scale!(similar(a), 5., a) == a*5 +@test scale!(similar(a), a, 5.) == a*5 +@test scale!(similar(a), a, [1.; 2.; 3.]) == a.*[1 2 3] +@test scale!(similar(a), [1.; 2.], a) == a.*[1; 2] +@test scale!(similar(a), a, [1; 2; 3]) == a.*[1 2 3] +@test scale!(similar(a), [1; 2], a) == a.*[1; 2] +@test_throws DimensionMismatch scale!(similar(a), a, ones(2)) +@test_throws DimensionMismatch scale!(similar(a), ones(3), a) +@test_throws DimensionMismatch scale!(Array(Float64, 3, 2), a, ones(3)) + +# scale real matrix by complex type +@test_throws InexactError scale!([1.0], 2.0im) +@test isequal(scale([1.0], 2.0im), Complex{Float64}[2.0im]) +@test isequal(scale(2.0im, [1.0]), Complex{Float64}[2.0im]) +@test isequal(scale(Float32[1.0], 2.0f0im), Complex{Float32}[2.0im]) +@test isequal(scale(Float32[1.0], 2.0im), Complex{Float64}[2.0im]) +@test isequal(scale(Float64[1.0], 2.0f0im), Complex{Float64}[2.0im]) +@test isequal(scale(Float32[1.0], big(2.0)im), Complex{BigFloat}[2.0im]) +@test isequal(scale(Float64[1.0], big(2.0)im), Complex{BigFloat}[2.0im]) +@test isequal(scale(BigFloat[1.0], 2.0im), Complex{BigFloat}[2.0im]) +@test isequal(scale(BigFloat[1.0], 2.0f0im), Complex{BigFloat}[2.0im]) diff --git a/test/linalg/lapack.jl b/test/linalg/lapack.jl index e272930a7d47f..f73c42325ad52 100644 --- a/test/linalg/lapack.jl +++ b/test/linalg/lapack.jl @@ -347,3 +347,28 @@ for elty in (Float32, Float64, Complex64, Complex128) X, rcond, f, b, r = LAPACK.gesvx!(C,D) @test_approx_eq X A\B end + +# Test our own linear algebra functionaly against LAPACK +for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) + for nn in (5,10,15) + if elty <: Real + A = convert(Matrix{elty}, randn(10,nn)) + else + A = convert(Matrix{elty}, complex(randn(10,nn),randn(10,nn))) + end ## LU (only equal for real because LAPACK uses different absolute value when choosing permutations) + if elty <: Real + FJulia = Base.LinAlg.generic_lufact!(copy(A)) + FLAPACK = Base.LinAlg.LAPACK.getrf!(copy(A)) + @test_approx_eq FJulia.factors FLAPACK[1] + @test_approx_eq FJulia.ipiv FLAPACK[2] + @test_approx_eq FJulia.info FLAPACK[3] + end + + ## QR + FJulia = invoke(qrfact!, Tuple{AbstractMatrix, Type{Val{false}}}, + copy(A), Val{false}) + FLAPACK = Base.LinAlg.LAPACK.geqrf!(copy(A)) + @test_approx_eq FJulia.factors FLAPACK[1] + @test_approx_eq FJulia.τ FLAPACK[2] + end +end diff --git a/test/linalg/matmul.jl b/test/linalg/matmul.jl new file mode 100644 index 0000000000000..4a2df859e7daf --- /dev/null +++ b/test/linalg/matmul.jl @@ -0,0 +1,134 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +using Base.Test + +## Test Julia fallbacks to BLAS routines + +# matrices with zero dimensions +for i = 1:10 + @test ones(0,5)*ones(5,3) == zeros(0,3) + @test ones(3,5)*ones(5,0) == zeros(3,0) + @test ones(3,0)*ones(0,4) == zeros(3,4) + @test ones(0,5)*ones(5,0) == zeros(0,0) + @test ones(0,0)*ones(0,4) == zeros(0,4) + @test ones(3,0)*ones(0,0) == zeros(3,0) + @test ones(0,0)*ones(0,0) == zeros(0,0) + @test Array(Float64, 5, 0) |> t -> t't == zeros(0,0) + @test Array(Float64, 5, 0) |> t -> t*t' == zeros(5,5) + @test Array(Complex128, 5, 0) |> t -> t't == zeros(0,0) + @test Array(Complex128, 5, 0) |> t -> t*t' == zeros(5,5) +end + +# 2x2 +A = [1 2; 3 4] +B = [5 6; 7 8] +@test A*B == [19 22; 43 50] +@test At_mul_B(A, B) == [26 30; 38 44] +@test A_mul_Bt(A, B) == [17 23; 39 53] +@test At_mul_Bt(A, B) == [23 31; 34 46] +Ai = A+(0.5*im).*B +Bi = B+(2.5*im).*A[[2,1],[2,1]] +@test Ai*Bi == [-21+53.5im -4.25+51.5im; -12+95.5im 13.75+85.5im] +@test Ac_mul_B(Ai, Bi) == [68.5-12im 57.5-28im; 88-3im 76.5-25im] +@test A_mul_Bc(Ai, Bi) == [64.5+5.5im 43+31.5im; 104-18.5im 80.5+31.5im] +@test Ac_mul_Bc(Ai, Bi) == [-28.25-66im 9.75-58im; -26-89im 21-73im] + +# 3x3 +A = [1 2 3; 4 5 6; 7 8 9].-5 +B = [1 0 5; 6 -10 3; 2 -4 -1] +@test A*B == [-26 38 -27; 1 -4 -6; 28 -46 15] +@test Ac_mul_B(A, B) == [-6 2 -25; 3 -12 -18; 12 -26 -11] +@test A_mul_Bc(A, B) == [-14 0 6; 4 -3 -3; 22 -6 -12] +@test Ac_mul_Bc(A, B) == [6 -8 -6; 12 -9 -9; 18 -10 -12] +Ai = A+(0.5*im).*B +Bi = B+(2.5*im).*A[[2,1,3],[2,3,1]] +@test Ai*Bi == [-44.75+13im 11.75-25im -38.25+30im; -47.75-16.5im -51.5+51.5im -56+6im; 16.75-4.5im -53.5+52im -15.5im] +@test Ac_mul_B(Ai, Bi) == [-21+2im -1.75+49im -51.25+19.5im; 25.5+56.5im -7-35.5im 22+35.5im; -3+12im -32.25+43im -34.75-2.5im] +@test A_mul_Bc(Ai, Bi) == [-20.25+15.5im -28.75-54.5im 22.25+68.5im; -12.25+13im -15.5+75im -23+27im; 18.25+im 1.5+94.5im -27-54.5im] +@test Ac_mul_Bc(Ai, Bi) == [1+2im 20.75+9im -44.75+42im; 19.5+17.5im -54-36.5im 51-14.5im; 13+7.5im 11.25+31.5im -43.25-14.5im] + +# Generic integer matrix multiplication +A = [1 2 3; 4 5 6] .- 3 +B = [2 -2; 3 -5; -4 7] +@test A*B == [-7 9; -4 9] +@test At_mul_Bt(A, B) == [-6 -11 15; -6 -13 18; -6 -15 21] +A = ones(Int, 2, 100) +B = ones(Int, 100, 3) +@test A*B == [100 100 100; 100 100 100] +A = rand(1:20, 5, 5) .- 10 +B = rand(1:20, 5, 5) .- 10 +@test At_mul_B(A, B) == A'*B +@test A_mul_Bt(A, B) == A*B' + +# Preallocated +C = Array(Int, size(A, 1), size(B, 2)) +@test A_mul_B!(C, A, B) == A*B +@test At_mul_B!(C, A, B) == A'*B +@test A_mul_Bt!(C, A, B) == A*B' +@test At_mul_Bt!(C, A, B) == A'*B' +v = [1,2,3] +C = Array(Int, 3, 3) +@test A_mul_Bt!(C, v, v) == v*v' +vf = map(Float64,v) +C = Array(Float64, 3, 3) +@test A_mul_Bt!(C, v, v) == v*v' + +# matrix algebra with subarrays of floats (stride != 1) +A = reshape(map(Float64,1:20),5,4) +Aref = A[1:2:end,1:2:end] +Asub = sub(A, 1:2:5, 1:2:4) +b = [1.2,-2.5] +@test (Aref*b) == (Asub*b) +@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref) +@test A_mul_Bt(Asub, Asub) == A_mul_Bt(Aref, Aref) +Ai = A .+ im +Aref = Ai[1:2:end,1:2:end] +Asub = sub(Ai, 1:2:5, 1:2:4) +@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref) +@test A_mul_Bc(Asub, Asub) == A_mul_Bc(Aref, Aref) + +# syrk & herk +A = reshape(1:1503, 501, 3).-750.0 +res = Float64[135228751 9979252 -115270247; 9979252 10481254 10983256; -115270247 10983256 137236759] +@test At_mul_B(A, A) == res +@test A_mul_Bt(A',A') == res +cutoff = 501 +A = reshape(1:6*cutoff,2*cutoff,3).-(6*cutoff)/2 +Asub = sub(A, 1:2:2*cutoff, 1:3) +Aref = A[1:2:2*cutoff, 1:3] +@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref) +Ai = A .- im +Asub = sub(Ai, 1:2:2*cutoff, 1:3) +Aref = Ai[1:2:2*cutoff, 1:3] +@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref) + +# matmul for types w/o sizeof (issue #1282) +A = Array(Complex{Int},10,10) +A[:] = complex(1,1) +A2 = A^2 +@test A2[1,1] == 20im + +# issue #6450 +@test dot(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 + +@test_throws DimensionMismatch dot([1.0,2.0,3.0], 1:2, [3.5,4.5,5.5], 1:3) +@test_throws BoundsError dot([1.0,2.0,3.0], 1:4, [3.5,4.5,5.5], 1:4) +@test_throws BoundsError dot([1.0,2.0,3.0], 1:3, [3.5,4.5,5.5], 2:4) +@test_throws DimensionMismatch dot(Complex128[1.0,2.0,3.0], 1:2, Complex128[3.5,4.5,5.5], 1:3) +@test_throws BoundsError dot(Complex128[1.0,2.0,3.0], 1:4, Complex128[3.5,4.5,5.5], 1:4) +@test_throws BoundsError dot(Complex128[1.0,2.0,3.0], 1:3, Complex128[3.5,4.5,5.5], 2:4) + +vecdot_(x,y) = invoke(vecdot, (Any,Any), x,y) # generic vecdot +let A = [1+2im 3+4im; 5+6im 7+8im], B = [2+7im 4+1im; 3+8im 6+5im] + @test vecdot(A,B) == dot(vec(A),vec(B)) == vecdot_(A,B) == vecdot(float(A),float(B)) + @test vecdot(Int[], Int[]) == 0 == vecdot_(Int[], Int[]) + @test_throws MethodError vecdot(Any[], Any[]) + @test_throws MethodError vecdot_(Any[], Any[]) + for n1 = 0:2, n2 = 0:2, d in (vecdot, vecdot_) + if n1 != n2 + @test_throws DimensionMismatch d(1:n1, 1:n2) + else + @test_approx_eq d(1:n1, 1:n2) vecnorm(1:n1)^2 + end + end +end diff --git a/test/linalg/qr.jl b/test/linalg/qr.jl new file mode 100644 index 0000000000000..a96b11967c3a6 --- /dev/null +++ b/test/linalg/qr.jl @@ -0,0 +1,126 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +debug = false +using Base.Test + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 +a2real = randn(n,n)/2 +a2img = randn(n,n)/2 +breal = randn(n,2)/2 +bimg = randn(n,2)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + + for eltyb in (Float32, Float64, Complex64, Complex128, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + +debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") +debug && println("QR decomposition (without pivoting)") + for i = 1:2 + let a = i == 1 ? a : sub(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : sub(b, 1:n - 1), n = i == 1 ? n : n - 1 + qra = @inferred qrfact(a) + @inferred qr(a) + q, r = qra[:Q], qra[:R] + @test_throws KeyError qra[:Z] + @test_approx_eq q'*full(q, thin = false) eye(n) + @test_approx_eq q*full(q, thin = false)' eye(n) + @test_approx_eq q*r a + @test_approx_eq_eps a*(qra\b) b 3000ε + @test_approx_eq full(qra) a + +debug && println("Thin QR decomposition (without pivoting)") + qra = @inferred qrfact(a[:,1:n1], Val{false}) + @inferred qr(a[:,1:n1], Val{false}) + q,r = qra[:Q], qra[:R] + @test_throws KeyError qra[:Z] + @test_approx_eq q'*full(q, thin=false) eye(n) + @test_approx_eq q'*full(q) eye(n, n1) + @test_approx_eq q*r a[:,1:n1] + @test_approx_eq_eps q*b[1:n1] full(q)*b[1:n1] 100ε + @test_approx_eq_eps q*b full(q, thin=false)*b 100ε + @test_throws DimensionMismatch q*b[1:n1 + 1] + +debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats + if eltya <: BlasFloat + @inferred qrfact(a, Val{true}) + @inferred qr(a, Val{true}) + end + qrpa = factorize(a[1:n1,:]) + q,r = qrpa[:Q], qrpa[:R] + @test_throws KeyError qrpa[:Z] + if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia + @test_approx_eq q'*full(q, thin=false) eye(n1) + @test_approx_eq q*full(q, thin=false)' eye(n1) + @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:n1,p] : a[1:n1,:] + @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:n1,:] + @test_approx_eq_eps a[1:n1,:]*(qrpa\b[1:n1]) b[1:n1] 5000ε + @test_approx_eq full(qrpa) a[1:5,:] + +debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats + qrpa = factorize(a[:,1:n1]) + q,r = qrpa[:Q], qrpa[:R] + @test_throws KeyError qrpa[:Z] + if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia + @test_approx_eq q'*full(q, thin=false) eye(n) + @test_approx_eq q*full(q, thin=false)' eye(n) + @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:n1] + @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:n1] + @test_approx_eq full(qrpa) a[:,1:5] + end + end + +debug && println("Matmul with QR factorizations") + if eltya != Int + qrpa = factorize(a[:,1:n1]) + q,r = qrpa[:Q], qrpa[:R] + @test_approx_eq A_mul_B!(full(q,thin=false)',q) eye(n) + @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q) + @test_approx_eq A_mul_Bc!(full(q,thin=false),q) eye(n) + @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q) + @test_throws BoundsError size(q,-1) + + qra = qrfact(a[:,1:n1], Val{false}) + q,r = qra[:Q], qra[:R] + @test_approx_eq A_mul_B!(full(q,thin=false)',q) eye(n) + @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q) + @test_approx_eq A_mul_Bc!(full(q,thin=false),q) eye(n) + @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q) + @test_throws BoundsError size(q,-1) + + @test_throws DimensionMismatch q * eye(Int8,n+4) + end + end +end + +# Because transpose(x) == x +@test_throws ErrorException transpose(qrfact(randn(3,3))) +@test_throws ErrorException ctranspose(qrfact(randn(3,3))) +@test_throws ErrorException transpose(qrfact(randn(3,3), Val{false})) +@test_throws ErrorException ctranspose(qrfact(randn(3,3), Val{false})) +@test_throws ErrorException transpose(qrfact(big(randn(3,3)))) +@test_throws ErrorException ctranspose(qrfact(big(randn(3,3)))) + +#Issue 7304 +let + A=[-√.5 -√.5; -√.5 √.5] + Q=full(qrfact(A)[:Q]) + @test vecnorm(A-Q) < eps() +end diff --git a/test/linalg/schur.jl b/test/linalg/schur.jl new file mode 100644 index 0000000000000..a632ce788b272 --- /dev/null +++ b/test/linalg/schur.jl @@ -0,0 +1,70 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +debug = false +using Base.Test + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + +debug && println("\ntype of a: ", eltya, "\n") +debug && println("Schur") + d,v = eig(a) + f = schurfact(a) + @test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a + @test_approx_eq sort(real(f[:values])) sort(real(d)) + @test_approx_eq sort(imag(f[:values])) sort(imag(d)) + @test istriu(f[:Schur]) || eltype(a)<:Real + @test_approx_eq full(f) a + @test_throws KeyError f[:A] + +debug && println("Reorder Schur") + # use asym for real schur to enforce tridiag structure + # avoiding partly selection of conj. eigenvalues + ordschura = eltya <: Complex ? a : asym + S = schurfact(ordschura) + select = bitrand(n) + O = ordschur(S, select) + sum(select) != 0 && @test_approx_eq S[:values][find(select)] O[:values][1:sum(select)] + @test_approx_eq O[:vectors]*O[:Schur]*O[:vectors]' ordschura + @test_throws KeyError f[:A] + +debug && println("Generalized Schur") + a1_sf = a[1:n1, 1:n1] + a2_sf = a[n1+1:n2, n1+1:n2] + f = schurfact(a1_sf, a2_sf) + @test_approx_eq f[:Q]*f[:S]*f[:Z]' a1_sf + @test_approx_eq f[:Q]*f[:T]*f[:Z]' a2_sf + @test istriu(f[:S]) || eltype(a)<:Real + @test istriu(f[:T]) || eltype(a)<:Real + @test_throws KeyError f[:A] + +debug && println("Reorder Generalized Schur") + a1_sf = a[1:n1, 1:n1] + a2_sf = a[n1+1:n2, n1+1:n2] + NS = schurfact(a1_sf, a2_sf) + # Currently just testing with selecting gen eig values < 1 + select = abs2(NS[:values]) .< 1 + m = sum(select) + S = ordschur(NS, select) + # Make sure that the new factorization stil factors matrix + @test_approx_eq S[:Q]*S[:S]*S[:Z]' a1_sf + @test_approx_eq S[:Q]*S[:T]*S[:Z]' a2_sf + # Make sure that we have sorted it correctly + @test_approx_eq NS[:values][find(select)] S[:values][1:m] +end diff --git a/test/linalg/special.jl b/test/linalg/special.jl new file mode 100644 index 0000000000000..af6c14564d6b5 --- /dev/null +++ b/test/linalg/special.jl @@ -0,0 +1,105 @@ +#This file is a part of Julia. License is MIT: http://julialang.org/license + +using Base.Test +debug = false + +n= 10 #Size of matrix to test +srand(1) + +debug && println("Test interconversion between special matrix types") +let a=[1.0:n;] + A=Diagonal(a) + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, LowerTriangular, UpperTriangular, Matrix] + debug && println("newtype is $(newtype)") + @test full(convert(newtype, A)) == full(A) + end + for newtype in [Base.LinAlg.UnitUpperTriangular, Base.LinAlg.UnitLowerTriangular] + @test_throws ArgumentError convert(newtype, A) + @test full(convert(newtype, Diagonal(ones(n)))) == eye(n) + end + + for isupper in (true, false) + debug && println("isupper is $(isupper)") + A=Bidiagonal(a, [1.0:n-1;], isupper) + for newtype in [Bidiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] + debug && println("newtype is $(newtype)") + @test full(convert(newtype, A)) == full(A) + @test full(newtype(A)) == full(A) + end + @test_throws ArgumentError convert(SymTridiagonal, A) + A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] + debug && println("newtype is $(newtype)") + @test full(convert(newtype, A)) == full(A) + @test full(newtype(A)) == full(A) + end + end + + A = SymTridiagonal(a, [1.0:n-1;]) + for newtype in [Tridiagonal, Matrix] + @test full(convert(newtype, A)) == full(A) + end + for newtype in [Diagonal, Bidiagonal] + @test_throws ArgumentError convert(newtype,A) + end + A = SymTridiagonal(a, zeros(n-1)) + @test full(convert(Bidiagonal,A)) == full(A) + + A = Tridiagonal(zeros(n-1), [1.0:n;], zeros(n-1)) #morally Diagonal + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix] + @test full(convert(newtype, A)) == full(A) + end + A = Tridiagonal(ones(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal + for newtype in [SymTridiagonal, Matrix] + @test full(convert(newtype, A)) == full(A) + end + for newtype in [Diagonal, Bidiagonal] + @test_throws ArgumentError convert(newtype,A) + end + A = Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal + @test full(convert(Bidiagonal, A)) == full(A) + A = UpperTriangular(Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1))) + @test full(convert(Bidiagonal, A)) == full(A) + A = Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)) #not morally Diagonal + @test full(convert(Bidiagonal, A)) == full(A) + A = LowerTriangular(Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1))) + @test full(convert(Bidiagonal, A)) == full(A) + + A = LowerTriangular(full(Diagonal(a))) #morally Diagonal + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix] + @test full(convert(newtype, A)) == full(A) + end + A = UpperTriangular(full(Diagonal(a))) #morally Diagonal + for newtype in [Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, Matrix] + @test full(convert(newtype, A)) == full(A) + end + A = UpperTriangular(triu(rand(n,n))) + for newtype in [Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal] + @test_throws ArgumentError convert(newtype,A) + end + A = Diagonal(a) + for newtype in [UpperTriangular, LowerTriangular] + @test full(convert(newtype,A)) == full(A) + end +end + +# Binary ops among special types +let a=[1.0:n;] + A=Diagonal(a) + Spectypes = [Diagonal, Bidiagonal, Tridiagonal, Matrix] + for (idx, type1) in enumerate(Spectypes) + for type2 in Spectypes + B = convert(type1,A) + C = convert(type2,A) + @test_approx_eq full(B + C) full(A + A) + @test_approx_eq full(B - C) full(A - A) + end + end + B = SymTridiagonal(a, ones(n-1)) + for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix] + @test_approx_eq full(B + convert(Spectype,A)) full(B + A) + @test_approx_eq full(convert(Spectype,A) + B) full(B + A) + @test_approx_eq full(B - convert(Spectype,A)) full(B - A) + @test_approx_eq full(convert(Spectype,A) - B) full(A - B) + end +end diff --git a/test/linalg/svd.jl b/test/linalg/svd.jl new file mode 100644 index 0000000000000..52bd0f1ea8f3b --- /dev/null +++ b/test/linalg/svd.jl @@ -0,0 +1,50 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +debug = false +using Base.Test + +using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted + +n = 10 + +# Split n into 2 parts for tests needing two matrices +n1 = div(n, 2) +n2 = 2*n1 + +srand(1234321) + +areal = randn(n,n)/2 +aimg = randn(n,n)/2 +a2real = randn(n,n)/2 +a2img = randn(n,n)/2 + +for eltya in (Float32, Float64, Complex64, Complex128, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + ε = εa = eps(abs(float(one(eltya)))) + +debug && println("\ntype of a: ", eltya, "\n") + +debug && println("singular value decomposition") + usv = svdfact(a) + @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a + @test_approx_eq full(usv) a + @test_approx_eq usv[:Vt]' usv[:V] + @test_throws KeyError usv[:Z] + b = rand(eltya,n) + @test_approx_eq usv\b a\b + +debug && println("Generalized svd") + a_svd = a[1:n1, :] + gsvd = svdfact(a,a_svd) + @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a + @test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a_svd + @test_approx_eq usv[:Vt]' usv[:V] + @test_throws KeyError usv[:Z] + α = eltya == Int ? -1 : rand(eltya) + β = svdfact(α) + @test β[:S] == [abs(α)] + @test svdvals(α) == [abs(α)] +end diff --git a/test/linalg/tridiag.jl b/test/linalg/tridiag.jl index 6f63ad8687b48..9a2081460b956 100644 --- a/test/linalg/tridiag.jl +++ b/test/linalg/tridiag.jl @@ -4,6 +4,120 @@ debug = false using Base.Test +# basic tridiagonal operations +n = 5 + +srand(123) + +d = 1 .+ rand(n) +dl = -rand(n-1) +du = -rand(n-1) +v = randn(n) +B = randn(n,2) + +for elty in (Float32, Float64, Complex64, Complex128, Int) + if elty == Int + srand(61516384) + d = rand(1:100, n) + dl = -rand(0:10, n-1) + du = -rand(0:10, n-1) + v = rand(1:100, n) + B = rand(1:100, n, 2) + else + d = convert(Vector{elty}, d) + dl = convert(Vector{elty}, dl) + du = convert(Vector{elty}, du) + v = convert(Vector{elty}, v) + B = convert(Matrix{elty}, B) + end + ε = eps(abs2(float(one(elty)))) + T = Tridiagonal(dl, d, du) + @test size(T, 1) == n + @test size(T) == (n, n) + F = diagm(d) + for i = 1:n-1 + F[i,i+1] = du[i] + F[i+1,i] = dl[i] + end + @test full(T) == F + + # elementary operations on tridiagonals + @test conj(T) == Tridiagonal(conj(dl), conj(d), conj(du)) + @test transpose(T) == Tridiagonal(du, d, dl) + @test ctranspose(T) == Tridiagonal(conj(du), conj(d), conj(dl)) + + # test interconversion of Tridiagonal and SymTridiagonal + @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl) + @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du) + @test SymTridiagonal(d, dl) + Tridiagonal(du, d, du) == SymTridiagonal(2d, dl+du) + + # tridiagonal linear algebra + @test_approx_eq T*v F*v + invFv = F\v + @test_approx_eq T\v invFv + # @test_approx_eq Base.solve(T,v) invFv + # @test_approx_eq Base.solve(T, B) F\B + Tlu = factorize(T) + x = Tlu\v + @test_approx_eq x invFv + @test_approx_eq det(T) det(F) + + @test_approx_eq T * Base.LinAlg.UnitUpperTriangular(eye(n)) F*eye(n) + @test_approx_eq T * Base.LinAlg.UnitLowerTriangular(eye(n)) F*eye(n) + @test_approx_eq T * UpperTriangular(eye(n)) F*eye(n) + @test_approx_eq T * LowerTriangular(eye(n)) F*eye(n) + + # symmetric tridiagonal + if elty <: Real + Ts = SymTridiagonal(d, dl) + Fs = full(Ts) + invFsv = Fs\v + Tldlt = ldltfact(Ts) + x = Tldlt\v + @test_approx_eq x invFsv + @test_approx_eq full(full(Tldlt)) Fs + end + + # eigenvalues/eigenvectors of symmetric tridiagonal + if elty === Float32 || elty === Float64 + DT, VT = eig(Ts) + D, Vecs = eig(Fs) + @test_approx_eq DT D + @test_approx_eq abs(VT'Vecs) eye(elty, n) + end + + # Test det(A::Matrix) + # In the long run, these tests should step through Strang's + # axiomatic definition of determinants. + # If all axioms are satisfied and all the composition rules work, + # all determinants will be correct except for floating point errors. + + # The determinant of the identity matrix should always be 1. + for i = 1:10 + A = eye(elty, i) + @test_approx_eq det(A) one(elty) + end + + # The determinant of a Householder reflection matrix should always be -1. + for i = 1:10 + A = eye(elty, 10) + A[i, i] = -one(elty) + @test_approx_eq det(A) -one(elty) + end + + # The determinant of a rotation matrix should always be 1. + if elty != Int + for theta = convert(Vector{elty}, pi ./ [1:4;]) + R = [cos(theta) -sin(theta); + sin(theta) cos(theta)] + @test_approx_eq convert(elty, det(R)) one(elty) + end + + # issue #1490 + @test_approx_eq_eps det(ones(elty, 3,3)) zero(elty) 3*eps(real(one(elty))) + end +end + #Test equivalence of eigenvectors/singular vectors taking into account possible phase (sign) differences function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S}, b::StridedVecOrMat{T}, error=nothing) n = size(a, 1) diff --git a/test/linalg1.jl b/test/linalg1.jl deleted file mode 100644 index e1ea207ad487a..0000000000000 --- a/test/linalg1.jl +++ /dev/null @@ -1,323 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -debug = false -using Base.Test - -using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted - -n = 10 - -# Split n into 2 parts for tests needing two matrices -n1 = div(n, 2) -n2 = 2*n1 - -srand(1234321) - -a = rand(n,n) -for elty in (Float32, Float64, Complex64, Complex128) - a = convert(Matrix{elty}, a) - # cond - @test_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01 - @test_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01 - @test_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01 - @test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01 -end - -areal = randn(n,n)/2 -aimg = randn(n,n)/2 -a2real = randn(n,n)/2 -a2img = randn(n,n)/2 -breal = randn(n,2)/2 -bimg = randn(n,2)/2 - -for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int) - a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) - a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) - asym = a'+a # symmetric indefinite - apd = a'*a # symmetric positive-definite - ε = εa = eps(abs(float(one(eltya)))) - - for eltyb in (Float32, Float64, Complex64, Complex128, Int) - b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal) - εb = eps(abs(float(one(eltyb)))) - ε = max(εa,εb) - -debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n") - -debug && println("(Automatic) Bunch-Kaufman factor of indefinite matrix") - if eltya != BigFloat && eltyb != BigFloat # Not implemented for BigFloat and I don't think it will. - bc1 = factorize(asym) - @test_approx_eq inv(bc1) * asym eye(n) - @test_approx_eq_eps asym * (bc1\b) b 1000ε - @test_approx_eq inv(a.' + a) * (a.' + a) eye(n) - @test size(bc1) == size(bc1.LD) - @test size(bc1,1) == size(bc1.LD,1) - @test size(bc1,2) == size(bc1.LD,2) - if eltya <: BlasReal - @test_throws ArgumentError bkfact(a) - end -debug && println("Bunch-Kaufman factors of a pos-def matrix") - bc2 = bkfact(apd) - @test_approx_eq inv(bc2) * apd eye(n) - @test_approx_eq_eps apd * (bc2\b) b 150000ε - @test ishermitian(bc2) == !issym(bc2) - - end - -debug && println("QR decomposition (without pivoting)") - for i = 1:2 - let a = i == 1 ? a : sub(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : sub(b, 1:n - 1), n = i == 1 ? n : n - 1 - qra = @inferred qrfact(a) - @inferred qr(a) - q, r = qra[:Q], qra[:R] - @test_throws KeyError qra[:Z] - @test_approx_eq q'*full(q, thin = false) eye(n) - @test_approx_eq q*full(q, thin = false)' eye(n) - @test_approx_eq q*r a - @test_approx_eq_eps a*(qra\b) b 3000ε - @test_approx_eq full(qra) a - -debug && println("Thin QR decomposition (without pivoting)") - qra = @inferred qrfact(a[:,1:n1], Val{false}) - @inferred qr(a[:,1:n1], Val{false}) - q,r = qra[:Q], qra[:R] - @test_throws KeyError qra[:Z] - @test_approx_eq q'*full(q, thin=false) eye(n) - @test_approx_eq q'*full(q) eye(n, n1) - @test_approx_eq q*r a[:,1:n1] - @test_approx_eq_eps q*b[1:n1] full(q)*b[1:n1] 100ε - @test_approx_eq_eps q*b full(q, thin=false)*b 100ε - @test_throws DimensionMismatch q*b[1:n1 + 1] - -debug && println("(Automatic) Fat (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats - if eltya <: BlasFloat - @inferred qrfact(a, Val{true}) - @inferred qr(a, Val{true}) - end - qrpa = factorize(a[1:n1,:]) - q,r = qrpa[:Q], qrpa[:R] - @test_throws KeyError qrpa[:Z] - if isa(qrpa,QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia - @test_approx_eq q'*full(q, thin=false) eye(n1) - @test_approx_eq q*full(q, thin=false)' eye(n1) - @test_approx_eq q*r isa(qrpa,QRPivoted) ? a[1:n1,p] : a[1:n1,:] - @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[1:n1,:] - @test_approx_eq_eps a[1:n1,:]*(qrpa\b[1:n1]) b[1:n1] 5000ε - @test_approx_eq full(qrpa) a[1:5,:] - -debug && println("(Automatic) Thin (pivoted) QR decomposition") # Pivoting is only implemented for BlasFloats - qrpa = factorize(a[:,1:n1]) - q,r = qrpa[:Q], qrpa[:R] - @test_throws KeyError qrpa[:Z] - if isa(qrpa, QRPivoted) p = qrpa[:p] end # Reconsider if pivoted QR gets implemented in julia - @test_approx_eq q'*full(q, thin=false) eye(n) - @test_approx_eq q*full(q, thin=false)' eye(n) - @test_approx_eq q*r isa(qrpa, QRPivoted) ? a[:,p] : a[:,1:n1] - @test_approx_eq isa(qrpa, QRPivoted) ? q*r[:,invperm(p)] : q*r a[:,1:n1] - @test_approx_eq full(qrpa) a[:,1:5] - end - end - -debug && println("Matmul with QR factorizations") - if eltya != Int - qrpa = factorize(a[:,1:n1]) - q,r = qrpa[:Q], qrpa[:R] - @test_approx_eq A_mul_B!(full(q,thin=false)',q) eye(n) - @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q) - @test_approx_eq A_mul_Bc!(full(q,thin=false),q) eye(n) - @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q) - @test_throws BoundsError size(q,-1) - - qra = qrfact(a[:,1:n1], Val{false}) - q,r = qra[:Q], qra[:R] - @test_approx_eq A_mul_B!(full(q,thin=false)',q) eye(n) - @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q) - @test_approx_eq A_mul_Bc!(full(q,thin=false),q) eye(n) - @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q) - @test_throws BoundsError size(q,-1) - - @test_throws DimensionMismatch q * eye(Int8,n+4) - end - -debug && println("non-symmetric eigen decomposition") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - d,v = eig(a) - for i in 1:size(a,2) @test_approx_eq a*v[:,i] d[i]*v[:,i] end - f = eigfact(a) - @test_approx_eq det(a) det(f) - @test_approx_eq inv(a) inv(f) - - num_fact = eigfact(one(eltya)) - @test num_fact.values[1] == one(eltya) - h = a + a' - @test_approx_eq minimum(eigvals(h)) eigmin(h) - @test_approx_eq maximum(eigvals(h)) eigmax(h) - @test_throws DomainError eigmin(a - a') - @test_throws DomainError eigmax(a - a') - end - -debug && println("symmetric generalized eigenproblem") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - asym_sg = asym[1:n1, 1:n1] - a_sg = a[:,n1+1:n2] - f = eigfact(asym_sg, a_sg'a_sg) - eig(asym_sg, a_sg'a_sg) # same result, but checks that method works - @test_approx_eq asym_sg*f[:vectors] scale(a_sg'a_sg*f[:vectors], f[:values]) - @test_approx_eq f[:values] eigvals(asym_sg, a_sg'a_sg) - @test_approx_eq_eps prod(f[:values]) prod(eigvals(asym_sg/(a_sg'a_sg))) 200ε - @test eigvecs(asym_sg, a_sg'a_sg) == f[:vectors] - @test_throws KeyError f[:Z] - end - -debug && println("Non-symmetric generalized eigenproblem") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - a1_nsg = a[1:n1, 1:n1] - a2_nsg = a[n1+1:n2, n1+1:n2] - f = eigfact(a1_nsg, a2_nsg) - eig(a1_nsg, a2_nsg) # same result, but checks that method works - @test_approx_eq a1_nsg*f[:vectors] scale(a2_nsg*f[:vectors], f[:values]) - @test_approx_eq f[:values] eigvals(a1_nsg, a2_nsg) - @test_approx_eq_eps prod(f[:values]) prod(eigvals(a1_nsg/a2_nsg)) 50000ε - @test eigvecs(a1_nsg, a2_nsg) == f[:vectors] - @test_throws KeyError f[:Z] - end - -debug && println("Schur") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - f = schurfact(a) - @test_approx_eq f[:vectors]*f[:Schur]*f[:vectors]' a - @test_approx_eq sort(real(f[:values])) sort(real(d)) - @test_approx_eq sort(imag(f[:values])) sort(imag(d)) - @test istriu(f[:Schur]) || eltype(a)<:Real - @test_approx_eq full(f) a - @test_throws KeyError f[:A] - end - -debug && println("Reorder Schur") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - # use asym for real schur to enforce tridiag structure - # avoiding partly selection of conj. eigenvalues - ordschura = eltya <: Complex ? a : asym - S = schurfact(ordschura) - select = bitrand(n) - O = ordschur(S, select) - sum(select) != 0 && @test_approx_eq S[:values][find(select)] O[:values][1:sum(select)] - @test_approx_eq O[:vectors]*O[:Schur]*O[:vectors]' ordschura - @test_throws KeyError f[:A] - end - -debug && println("Generalized Schur") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - a1_sf = a[1:n1, 1:n1] - a2_sf = a[n1+1:n2, n1+1:n2] - f = schurfact(a1_sf, a2_sf) - @test_approx_eq f[:Q]*f[:S]*f[:Z]' a1_sf - @test_approx_eq f[:Q]*f[:T]*f[:Z]' a2_sf - @test istriu(f[:S]) || eltype(a)<:Real - @test istriu(f[:T]) || eltype(a)<:Real - @test_throws KeyError f[:A] - end - -debug && println("Reorder Generalized Schur") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in Julia - a1_sf = a[1:n1, 1:n1] - a2_sf = a[n1+1:n2, n1+1:n2] - NS = schurfact(a1_sf, a2_sf) - # Currently just testing with selecting gen eig values < 1 - select = abs2(NS[:values]) .< 1 - m = sum(select) - S = ordschur(NS, select) - # Make sure that the new factorization stil factors matrix - @test_approx_eq S[:Q]*S[:S]*S[:Z]' a1_sf - @test_approx_eq S[:Q]*S[:T]*S[:Z]' a2_sf - # Make sure that we have sorted it correctly - @test_approx_eq NS[:values][find(select)] S[:values][1:m] - end - -debug && println("singular value decomposition") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - usv = svdfact(a) - @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a - @test_approx_eq full(usv) a - @test_approx_eq usv[:Vt]' usv[:V] - @test_throws KeyError usv[:Z] - b = rand(eltya,n) - @test_approx_eq usv\b a\b - end - -debug && println("Generalized svd") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - a_svd = a[1:n1, :] - gsvd = svdfact(a,a_svd) - @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a - @test_approx_eq gsvd[:V]*gsvd[:D2]*gsvd[:R]*gsvd[:Q]' a_svd - @test_approx_eq usv[:Vt]' usv[:V] - @test_throws KeyError usv[:Z] - α = eltya == Int ? -1 : rand(eltya) - β = svdfact(α) - @test β[:S] == [abs(α)] - @test svdvals(α) == [abs(α)] - end - -debug && println("Solve square general system of equations") - κ = cond(a,1) - x = a \ b - @test_throws DimensionMismatch b'\b - @test_throws DimensionMismatch b\b' - @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit! - -debug && println("Test nullspace") - if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia - a15null = nullspace(a[:,1:n1]') - @test rank([a[:,1:n1] a15null]) == 10 - @test_approx_eq_eps norm(a[:,1:n1]'a15null, Inf) zero(eltya) 300ε - @test_approx_eq_eps norm(a15null'a[:,1:n1], Inf) zero(eltya) 400ε - @test size(nullspace(b), 2) == 0 - end - end # for eltyb - -debug && println("\ntype of a: ", eltya, "\n") - -debug && println("Test pinv") - if eltya != BigFloat # Revisit when implemented in julia - pinva15 = pinv(a[:,1:n1]) - @test_approx_eq a[:,1:n1]*pinva15*a[:,1:n1] a[:,1:n1] - @test_approx_eq pinva15*a[:,1:n1]*pinva15 pinva15 - end - - # if isreal(a) -debug && println("Matrix square root") - if eltya != BigFloat # Revisit when implemented in julia - asq = sqrtm(a) - @test_approx_eq asq*asq a - asymsq = sqrtm(asym) - @test_approx_eq asymsq*asymsq asym - end - -debug && println("Lyapunov/Sylvester") - if eltya != BigFloat - let - x = lyap(a, a2) - @test_approx_eq -a2 a*x + x*a' - x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n]) - @test_approx_eq -a2[1:3, 4:n] a[1:3, 1:3]*x2 + x2*a[4:n, 4:n] - end - end -end # for eltya - -#6941 -#@test (ones(10^7,4)*ones(4))[3] == 4.0 - -# test diff, throw ArgumentError for invalid dimension argument -let X = [3 9 5; - 7 4 2; - 2 1 10] - @test diff(X,1) == [4 -5 -3; -5 -3 8] - @test diff(X,2) == [6 -4; -3 -2; -1 9] - @test_throws ArgumentError diff(X,3) - @test_throws ArgumentError diff(X,-1) -end - -x = float([1:12;]) -y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99] -@test_approx_eq linreg(x,y) [2.5559090909090867, 1.6960139860139862] diff --git a/test/linalg2.jl b/test/linalg2.jl deleted file mode 100644 index 4bdfeef1d588c..0000000000000 --- a/test/linalg2.jl +++ /dev/null @@ -1,395 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -debug = false - -import Base.LinAlg -import Base.LinAlg: BlasComplex, BlasFloat, BlasReal, UnitUpperTriangular, UnitLowerTriangular - -# basic tridiagonal operations -n = 5 - -srand(123) - -d = 1 .+ rand(n) -dl = -rand(n-1) -du = -rand(n-1) -v = randn(n) -B = randn(n,2) - -for elty in (Float32, Float64, Complex64, Complex128, Int) - if elty == Int - srand(61516384) - d = rand(1:100, n) - dl = -rand(0:10, n-1) - du = -rand(0:10, n-1) - v = rand(1:100, n) - B = rand(1:100, n, 2) - else - d = convert(Vector{elty}, d) - dl = convert(Vector{elty}, dl) - du = convert(Vector{elty}, du) - v = convert(Vector{elty}, v) - B = convert(Matrix{elty}, B) - end - ε = eps(abs2(float(one(elty)))) - T = Tridiagonal(dl, d, du) - @test size(T, 1) == n - @test size(T) == (n, n) - F = diagm(d) - for i = 1:n-1 - F[i,i+1] = du[i] - F[i+1,i] = dl[i] - end - @test full(T) == F - - # elementary operations on tridiagonals - @test conj(T) == Tridiagonal(conj(dl), conj(d), conj(du)) - @test transpose(T) == Tridiagonal(du, d, dl) - @test ctranspose(T) == Tridiagonal(conj(du), conj(d), conj(dl)) - - # test interconversion of Tridiagonal and SymTridiagonal - @test Tridiagonal(dl, d, dl) == SymTridiagonal(d, dl) - @test Tridiagonal(dl, d, du) + Tridiagonal(du, d, dl) == SymTridiagonal(2d, dl+du) - @test SymTridiagonal(d, dl) + Tridiagonal(du, d, du) == SymTridiagonal(2d, dl+du) - - # tridiagonal linear algebra - @test_approx_eq T*v F*v - invFv = F\v - @test_approx_eq T\v invFv - # @test_approx_eq Base.solve(T,v) invFv - # @test_approx_eq Base.solve(T, B) F\B - Tlu = factorize(T) - x = Tlu\v - @test_approx_eq x invFv - @test_approx_eq det(T) det(F) - - @test_approx_eq T * UnitUpperTriangular(eye(n)) F*eye(n) - @test_approx_eq T * UnitLowerTriangular(eye(n)) F*eye(n) - @test_approx_eq T * UpperTriangular(eye(n)) F*eye(n) - @test_approx_eq T * LowerTriangular(eye(n)) F*eye(n) - - # symmetric tridiagonal - if elty <: Real - Ts = SymTridiagonal(d, dl) - Fs = full(Ts) - invFsv = Fs\v - Tldlt = ldltfact(Ts) - x = Tldlt\v - @test_approx_eq x invFsv - @test_approx_eq full(full(Tldlt)) Fs - end - - # eigenvalues/eigenvectors of symmetric tridiagonal - if elty === Float32 || elty === Float64 - DT, VT = eig(Ts) - D, Vecs = eig(Fs) - @test_approx_eq DT D - @test_approx_eq abs(VT'Vecs) eye(elty, n) - end - - # Test det(A::Matrix) - # In the long run, these tests should step through Strang's - # axiomatic definition of determinants. - # If all axioms are satisfied and all the composition rules work, - # all determinants will be correct except for floating point errors. - - # The determinant of the identity matrix should always be 1. - for i = 1:10 - A = eye(elty, i) - @test_approx_eq det(A) one(elty) - end - - # The determinant of a Householder reflection matrix should always be -1. - for i = 1:10 - A = eye(elty, 10) - A[i, i] = -one(elty) - @test_approx_eq det(A) -one(elty) - end - - # The determinant of a rotation matrix should always be 1. - if elty != Int - for theta = convert(Vector{elty}, pi ./ [1:4;]) - R = [cos(theta) -sin(theta); - sin(theta) cos(theta)] - @test_approx_eq convert(elty, det(R)) one(elty) - end - - # issue #1490 - @test_approx_eq_eps det(ones(elty, 3,3)) zero(elty) 3*eps(real(one(elty))) - end -end - -# Generic BLAS tests -srand(100) -# syr2k! and her2k! -for elty in (Float32, Float64, Complex64, Complex128) - U = randn(5,2) - V = randn(5,2) - if elty == Complex64 || elty == Complex128 - U = complex(U, U) - V = complex(V, V) - end - U = convert(Array{elty, 2}, U) - V = convert(Array{elty, 2}, V) - @test_approx_eq tril(LinAlg.BLAS.syr2k('L','N',U,V)) tril(U*V.' + V*U.') - @test_approx_eq triu(LinAlg.BLAS.syr2k('U','N',U,V)) triu(U*V.' + V*U.') - @test_approx_eq tril(LinAlg.BLAS.syr2k('L','T',U,V)) tril(U.'*V + V.'*U) - @test_approx_eq triu(LinAlg.BLAS.syr2k('U','T',U,V)) triu(U.'*V + V.'*U) -end - -for elty in (Complex64, Complex128) - U = randn(5,2) - V = randn(5,2) - if elty == Complex64 || elty == Complex128 - U = complex(U, U) - V = complex(V, V) - end - U = convert(Array{elty, 2}, U) - V = convert(Array{elty, 2}, V) - @test_approx_eq tril(LinAlg.BLAS.her2k('L','N',U,V)) tril(U*V' + V*U') - @test_approx_eq triu(LinAlg.BLAS.her2k('U','N',U,V)) triu(U*V' + V*U') - @test_approx_eq tril(LinAlg.BLAS.her2k('L','C',U,V)) tril(U'*V + V'*U) - @test_approx_eq triu(LinAlg.BLAS.her2k('U','C',U,V)) triu(U'*V + V'*U) -end - -# Test gradient -for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128) - if elty <: Real - x = convert(Vector{elty}, [1:3;]) - g = ones(elty, 3) - else - x = convert(Vector{elty}, complex([1:3;], [1:3;])) - g = convert(Vector{elty}, complex(ones(3), ones(3))) - end - @test_approx_eq gradient(x) g -end - -# Test our own linear algebra functionaly against LAPACK -for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) - for nn in (5,10,15) - if elty <: Real - A = convert(Matrix{elty}, randn(10,nn)) - else - A = convert(Matrix{elty}, complex(randn(10,nn),randn(10,nn))) - end ## LU (only equal for real because LAPACK uses different absolute value when choosing permutations) - if elty <: Real - FJulia = Base.LinAlg.generic_lufact!(copy(A)) - FLAPACK = Base.LinAlg.LAPACK.getrf!(copy(A)) - @test_approx_eq FJulia.factors FLAPACK[1] - @test_approx_eq FJulia.ipiv FLAPACK[2] - @test_approx_eq FJulia.info FLAPACK[3] - end - - ## QR - FJulia = invoke(qrfact!, Tuple{AbstractMatrix, Type{Val{false}}}, - copy(A), Val{false}) - FLAPACK = Base.LinAlg.LAPACK.geqrf!(copy(A)) - @test_approx_eq FJulia.factors FLAPACK[1] - @test_approx_eq FJulia.τ FLAPACK[2] - end -end - -# Tests norms -nnorm = 10 -mmat = 10 -nmat = 8 -for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt) - debug && println(elty) - - ## Vector - x = ones(elty,10) - xs = sub(x,1:2:10) - @test_approx_eq norm(x, -Inf) 1 - @test_approx_eq norm(x, -1) 1/10 - @test_approx_eq norm(x, 0) 10 - @test_approx_eq norm(x, 1) 10 - @test_approx_eq norm(x, 2) sqrt(10) - @test_approx_eq norm(x, 3) cbrt(10) - @test_approx_eq norm(x, Inf) 1 - @test_approx_eq norm(xs, -Inf) 1 - @test_approx_eq norm(xs, -1) 1/5 - @test_approx_eq norm(xs, 0) 5 - @test_approx_eq norm(xs, 1) 5 - @test_approx_eq norm(xs, 2) sqrt(5) - @test_approx_eq norm(xs, 3) cbrt(5) - @test_approx_eq norm(xs, Inf) 1 - - ## Number - norm(x[1:1]) === norm(x[1], -Inf) - norm(x[1:1]) === norm(x[1], 0) - norm(x[1:1]) === norm(x[1], 1) - norm(x[1:1]) === norm(x[1], 2) - norm(x[1:1]) === norm(x[1], Inf) - - for i = 1:10 - x = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : - elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : - convert(Vector{elty}, randn(nnorm)) - xs = sub(x,1:2:nnorm) - y = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) : - elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) : - convert(Vector{elty}, randn(nnorm)) - ys = sub(y,1:2:nnorm) - α = elty <: Integer ? randn() : - elty <: Complex ? convert(elty, complex(randn(),randn())) : - convert(elty, randn()) - # Absolute homogeneity - @test_approx_eq norm(α*x,-Inf) abs(α)*norm(x,-Inf) - @test_approx_eq norm(α*x,-1) abs(α)*norm(x,-1) - @test_approx_eq norm(α*x,1) abs(α)*norm(x,1) - @test_approx_eq norm(α*x) abs(α)*norm(x) # two is default - @test_approx_eq norm(α*x,3) abs(α)*norm(x,3) - @test_approx_eq norm(α*x,Inf) abs(α)*norm(x,Inf) - - @test_approx_eq norm(α*xs,-Inf) abs(α)*norm(xs,-Inf) - @test_approx_eq norm(α*xs,-1) abs(α)*norm(xs,-1) - @test_approx_eq norm(α*xs,1) abs(α)*norm(xs,1) - @test_approx_eq norm(α*xs) abs(α)*norm(xs) # two is default - @test_approx_eq norm(α*xs,3) abs(α)*norm(xs,3) - @test_approx_eq norm(α*xs,Inf) abs(α)*norm(xs,Inf) - - # Triangle inequality - @test norm(x + y,1) <= norm(x,1) + norm(y,1) - @test norm(x + y) <= norm(x) + norm(y) # two is default - @test norm(x + y,3) <= norm(x,3) + norm(y,3) - @test norm(x + y,Inf) <= norm(x,Inf) + norm(y,Inf) - - @test norm(xs + ys,1) <= norm(xs,1) + norm(ys,1) - @test norm(xs + ys) <= norm(xs) + norm(ys) # two is default - @test norm(xs + ys,3) <= norm(xs,3) + norm(ys,3) - @test norm(xs + ys,Inf) <= norm(xs,Inf) + norm(ys,Inf) - - # Against vectorized versions - @test_approx_eq norm(x,-Inf) minimum(abs(x)) - @test_approx_eq norm(x,-1) inv(sum(1./abs(x))) - @test_approx_eq norm(x,0) sum(x .!= 0) - @test_approx_eq norm(x,1) sum(abs(x)) - @test_approx_eq norm(x) sqrt(sum(abs2(x))) - @test_approx_eq norm(x,3) cbrt(sum(abs(x).^3.)) - @test_approx_eq norm(x,Inf) maximum(abs(x)) - end - ## Matrix (Operator) - A = ones(elty,10,10) - As = sub(A,1:5,1:5) - @test_approx_eq norm(A, 1) 10 - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(A, 2) 10 - @test_approx_eq norm(A, Inf) 10 - @test_approx_eq norm(As, 1) 5 - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(As, 2) 5 - @test_approx_eq norm(As, Inf) 5 - - for i = 1:10 - A = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : - elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : - convert(Matrix{elty}, randn(mmat, nmat)) - As = sub(A,1:nmat,1:nmat) - B = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) : - elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) : - convert(Matrix{elty}, randn(mmat, nmat)) - Bs = sub(B,1:nmat,1:nmat) - α = elty <: Integer ? randn() : - elty <: Complex ? convert(elty, complex(randn(),randn())) : - convert(elty, randn()) - - # Absolute homogeneity - @test_approx_eq norm(α*A,1) abs(α)*norm(A,1) - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*A) abs(α)*norm(A) # two is default - @test_approx_eq norm(α*A,Inf) abs(α)*norm(A,Inf) - - @test_approx_eq norm(α*As,1) abs(α)*norm(As,1) - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*As) abs(α)*norm(As) # two is default - @test_approx_eq norm(α*As,Inf) abs(α)*norm(As,Inf) - - # Triangle inequality - @test norm(A + B,1) <= norm(A,1) + norm(B,1) - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A + B) <= norm(A) + norm(B) # two is default - @test norm(A + B,Inf) <= norm(A,Inf) + norm(B,Inf) - - @test norm(As + Bs,1) <= norm(As,1) + norm(Bs,1) - elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(As + Bs) <= norm(As) + norm(Bs) # two is default - @test norm(As + Bs,Inf) <= norm(As,Inf) + norm(Bs,Inf) - - # vecnorm: - for p = -2:3 - @test norm(reshape(A, length(A)), p) == vecnorm(A, p) - end - - # issue #10234 - if elty <: AbstractFloat || elty <: Complex - let z = zeros(elty, 100) - z[1] = -Inf - for p in [-2,-1.5,-1,-0.5,0.5,1,1.5,2,Inf] - @test norm(z, p) == (p < 0 ? 0 : Inf) - @test norm(elty[Inf],p) == Inf - end - end - end - end -end - -# issue #10234 -@test norm(Any[Inf],-2) == norm(Any[Inf],-1) == norm(Any[Inf],1) == norm(Any[Inf],1.5) == norm(Any[Inf],2) == norm(Any[Inf],Inf) == Inf - -# overflow/underflow in norms: -@test_approx_eq norm(Float64[1e-300, 1], -3)*1e300 1 -@test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1 - -## Issue related tests -# issue #1447 -let - A = [1.+0.im 0; 0 1] - B = pinv(A) - for i = 1:4 - @test_approx_eq A[i] B[i] - end -end - -# issue #2246 -let - A = [1 2 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0] - Asq = sqrtm(A) - @test_approx_eq Asq*Asq A - A2 = sub(A, 1:2, 1:2) - A2sq = sqrtm(A2) - @test_approx_eq A2sq*A2sq A2 -end - -let - N = 3 - @test_approx_eq log(det(eye(N))) logdet(eye(N)) -end - -# issue #2637 -let - a = [1, 2, 3] - b = [4, 5, 6] - @test kron(eye(2),eye(2)) == eye(4) - @test kron(a,b) == [4,5,6,8,10,12,12,15,18] - @test kron(a',b') == [4 5 6 8 10 12 12 15 18] - @test kron(a,b') == [4 5 6; 8 10 12; 12 15 18] - @test kron(a',b) == [4 8 12; 5 10 15; 6 12 18] - @test kron(a,eye(2)) == [1 0; 0 1; 2 0; 0 2; 3 0; 0 3] - @test kron(eye(2),a) == [ 1 0; 2 0; 3 0; 0 1; 0 2; 0 3] - @test kron(eye(2),2) == 2*eye(2) - @test kron(3,eye(3)) == 3*eye(3) - @test kron(a,2) == [2, 4, 6] - @test kron(b',2) == [8 10 12] -end - -# issue #4796 -let - dim=2 - S=zeros(Complex,dim,dim) - T=zeros(Complex,dim,dim) - T[:] = 1 - z = 2.5 + 1.5im - S[1] = z - @test S*T == [z z; 0 0] -end - -#Issue 7304 -let - A=[-√.5 -√.5; -√.5 √.5] - Q=full(qrfact(A)[:Q]) - @test vecnorm(A-Q) < eps() -end diff --git a/test/linalg3.jl b/test/linalg3.jl deleted file mode 100644 index 51f5605da2941..0000000000000 --- a/test/linalg3.jl +++ /dev/null @@ -1,333 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -using Base.Test - -## Least squares solutions -a = [ones(20) 1:20 1:20] -b = reshape(eye(8, 5), 20, 2) -for elty in (Float32, Float64, Complex64, Complex128) - a = convert(Matrix{elty}, a) - b = convert(Matrix{elty}, b) - - # Vector rhs - x = a[:,1:2]\b[:,1] - @test_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384) - - # Matrix rhs - x = a[:,1:2]\b - @test_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031) - - # Rank deficient - x = a\b - @test_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031) - - # Underdetermined minimum norm - x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1]) - @test_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5]) - - # symmetric, positive definite - @test_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3]) - - # symmetric, indefinite - @test_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3) -end - -## Test Julia fallbacks to BLAS routines - -# matrices with zero dimensions -for i = 1:10 - @test ones(0,5)*ones(5,3) == zeros(0,3) - @test ones(3,5)*ones(5,0) == zeros(3,0) - @test ones(3,0)*ones(0,4) == zeros(3,4) - @test ones(0,5)*ones(5,0) == zeros(0,0) - @test ones(0,0)*ones(0,4) == zeros(0,4) - @test ones(3,0)*ones(0,0) == zeros(3,0) - @test ones(0,0)*ones(0,0) == zeros(0,0) - @test Array(Float64, 5, 0) |> t -> t't == zeros(0,0) - @test Array(Float64, 5, 0) |> t -> t*t' == zeros(5,5) - @test Array(Complex128, 5, 0) |> t -> t't == zeros(0,0) - @test Array(Complex128, 5, 0) |> t -> t*t' == zeros(5,5) -end - -# 2x2 -A = [1 2; 3 4] -B = [5 6; 7 8] -@test A*B == [19 22; 43 50] -@test At_mul_B(A, B) == [26 30; 38 44] -@test A_mul_Bt(A, B) == [17 23; 39 53] -@test At_mul_Bt(A, B) == [23 31; 34 46] -Ai = A+(0.5*im).*B -Bi = B+(2.5*im).*A[[2,1],[2,1]] -@test Ai*Bi == [-21+53.5im -4.25+51.5im; -12+95.5im 13.75+85.5im] -@test Ac_mul_B(Ai, Bi) == [68.5-12im 57.5-28im; 88-3im 76.5-25im] -@test A_mul_Bc(Ai, Bi) == [64.5+5.5im 43+31.5im; 104-18.5im 80.5+31.5im] -@test Ac_mul_Bc(Ai, Bi) == [-28.25-66im 9.75-58im; -26-89im 21-73im] - -# 3x3 -A = [1 2 3; 4 5 6; 7 8 9].-5 -B = [1 0 5; 6 -10 3; 2 -4 -1] -@test A*B == [-26 38 -27; 1 -4 -6; 28 -46 15] -@test Ac_mul_B(A, B) == [-6 2 -25; 3 -12 -18; 12 -26 -11] -@test A_mul_Bc(A, B) == [-14 0 6; 4 -3 -3; 22 -6 -12] -@test Ac_mul_Bc(A, B) == [6 -8 -6; 12 -9 -9; 18 -10 -12] -Ai = A+(0.5*im).*B -Bi = B+(2.5*im).*A[[2,1,3],[2,3,1]] -@test Ai*Bi == [-44.75+13im 11.75-25im -38.25+30im; -47.75-16.5im -51.5+51.5im -56+6im; 16.75-4.5im -53.5+52im -15.5im] -@test Ac_mul_B(Ai, Bi) == [-21+2im -1.75+49im -51.25+19.5im; 25.5+56.5im -7-35.5im 22+35.5im; -3+12im -32.25+43im -34.75-2.5im] -@test A_mul_Bc(Ai, Bi) == [-20.25+15.5im -28.75-54.5im 22.25+68.5im; -12.25+13im -15.5+75im -23+27im; 18.25+im 1.5+94.5im -27-54.5im] -@test Ac_mul_Bc(Ai, Bi) == [1+2im 20.75+9im -44.75+42im; 19.5+17.5im -54-36.5im 51-14.5im; 13+7.5im 11.25+31.5im -43.25-14.5im] - -# Generic integer matrix multiplication -A = [1 2 3; 4 5 6] .- 3 -B = [2 -2; 3 -5; -4 7] -@test A*B == [-7 9; -4 9] -@test At_mul_Bt(A, B) == [-6 -11 15; -6 -13 18; -6 -15 21] -A = ones(Int, 2, 100) -B = ones(Int, 100, 3) -@test A*B == [100 100 100; 100 100 100] -A = rand(1:20, 5, 5) .- 10 -B = rand(1:20, 5, 5) .- 10 -@test At_mul_B(A, B) == A'*B -@test A_mul_Bt(A, B) == A*B' - -# Preallocated -C = Array(Int, size(A, 1), size(B, 2)) -@test A_mul_B!(C, A, B) == A*B -@test At_mul_B!(C, A, B) == A'*B -@test A_mul_Bt!(C, A, B) == A*B' -@test At_mul_Bt!(C, A, B) == A'*B' -v = [1,2,3] -C = Array(Int, 3, 3) -@test A_mul_Bt!(C, v, v) == v*v' -vf = map(Float64,v) -C = Array(Float64, 3, 3) -@test A_mul_Bt!(C, v, v) == v*v' - -# matrix algebra with subarrays of floats (stride != 1) -A = reshape(map(Float64,1:20),5,4) -Aref = A[1:2:end,1:2:end] -Asub = sub(A, 1:2:5, 1:2:4) -b = [1.2,-2.5] -@test (Aref*b) == (Asub*b) -@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref) -@test A_mul_Bt(Asub, Asub) == A_mul_Bt(Aref, Aref) -Ai = A .+ im -Aref = Ai[1:2:end,1:2:end] -Asub = sub(Ai, 1:2:5, 1:2:4) -@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref) -@test A_mul_Bc(Asub, Asub) == A_mul_Bc(Aref, Aref) - -# syrk & herk -A = reshape(1:1503, 501, 3).-750.0 -res = Float64[135228751 9979252 -115270247; 9979252 10481254 10983256; -115270247 10983256 137236759] -@test At_mul_B(A, A) == res -@test A_mul_Bt(A',A') == res -cutoff = 501 -A = reshape(1:6*cutoff,2*cutoff,3).-(6*cutoff)/2 -Asub = sub(A, 1:2:2*cutoff, 1:3) -Aref = A[1:2:2*cutoff, 1:3] -@test At_mul_B(Asub, Asub) == At_mul_B(Aref, Aref) -Ai = A .- im -Asub = sub(Ai, 1:2:2*cutoff, 1:3) -Aref = Ai[1:2:2*cutoff, 1:3] -@test Ac_mul_B(Asub, Asub) == Ac_mul_B(Aref, Aref) - -# Matrix exponential -for elty in (Float32, Float64, Complex64, Complex128) - A1 = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4]) - eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181 127.781085523182; - 183.765138646367 183.765138646366 163.679601723179; - 71.797032399996 91.8825693231832 111.968106246371]') - @test_approx_eq expm(A1) eA1 - - A2 = convert(Matrix{elty}, - [29.87942128909879 0.7815750847907159 -2.289519314033932; - 0.7815750847907159 25.72656945571064 8.680737820540137; - -2.289519314033932 8.680737820540137 34.39400925519054]) - eA2 = convert(Matrix{elty}, - [ 5496313853692458.0 -18231880972009236.0 -30475770808580460.0; - -18231880972009252.0 60605228702221920.0 101291842930249760.0; - -30475770808580480.0 101291842930249728.0 169294411240851968.0]) - @test_approx_eq expm(A2) eA2 - - A3 = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52]) - eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812 -4.934938326092; - 0.367879439109187 1.47151775849686 1.10363831732856; - 0.135335281175235 0.406005843524598 0.541341126763207]') - @test_approx_eq expm(A3) eA3 - - # Hessenberg - @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty}, - [4.000000000000000 -1.414213562373094 -1.414213562373095 - -1.414213562373095 4.999999999999996 -0.000000000000000 - 0 -0.000000000000002 3.000000000000000]) - @test_throws KeyError hessfact(A1)[:Z] - @test_approx_eq full(hessfact(A1)) A1 - @test_approx_eq full(Base.LinAlg.HessenbergQ(hessfact(A1))) full(hessfact(A1)[:Q]) -end - -for elty in (Float64, Complex{Float64}) - - A4 = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps(); - 1/3 1/4 1/5 1/6; - 1/4 1/5 1/6 1/7; - 1/5 1/6 1/7 1/8]) - @test_approx_eq expm(logm(A4)) A4 - - OLD_STDERR = STDERR - rd,wr = redirect_stderr() - A5 = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1]) - @test_approx_eq expm(logm(A5)) A5 - s = readline(rd) - @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.") - - A6 = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11]) - @test_approx_eq expm(logm(A6)) A6 - s = readline(rd) - @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprimary matrix logarithm will be returned.") - redirect_stderr(OLD_STDERR) - - A7 = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]) - @test_approx_eq expm(logm(A7)) A7 - -end - -A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1]; -@test_approx_eq expm(logm(A8)) A8 - -# issue 5116 -A9 = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0] -eA9 = [-0.999786072879326 -0.065407069689389 0.0 0.0 - 0.006540706968939 -0.999786072879326 0.0 0.0 - 0.0 0.0 1.0 0.0 - 0.013081413937878 -3.999572145758650 0.0 1.0] -@test_approx_eq expm(A9) eA9 - -# issue 5116 -A10 = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.] -eA10 = [ 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im - 0.0+0.0im 1.543080634815244+0.0im 0.0-1.175201193643801im 0.0+0.0im - 0.0+0.0im 0.0+1.175201193643801im 1.543080634815243+0.0im 0.0+0.0im - 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im] -@test_approx_eq expm(A10) eA10 - -# matmul for types w/o sizeof (issue #1282) -A = Array(Complex{Int},10,10) -A[:] = complex(1,1) -A2 = A^2 -@test A2[1,1] == 20im - -# test sparse matrix norms -Ac = sprandn(10,10,.1) + im* sprandn(10,10,.1) -Ar = sprandn(10,10,.1) -Ai = ceil(Int,Ar*100) -@test_approx_eq norm(Ac,1) norm(full(Ac),1) -@test_approx_eq norm(Ac,Inf) norm(full(Ac),Inf) -@test_approx_eq vecnorm(Ac) vecnorm(full(Ac)) -@test_approx_eq norm(Ar,1) norm(full(Ar),1) -@test_approx_eq norm(Ar,Inf) norm(full(Ar),Inf) -@test_approx_eq vecnorm(Ar) vecnorm(full(Ar)) -@test_approx_eq norm(Ai,1) norm(full(Ai),1) -@test_approx_eq norm(Ai,Inf) norm(full(Ai),Inf) -@test_approx_eq vecnorm(Ai) vecnorm(full(Ai)) - -# 2-argument version of scale -a = reshape([1.:6;], (2,3)) -@test scale(a, 5.) == a*5 -@test scale(5., a) == a*5 -@test scale(a, [1.; 2.; 3.]) == a.*[1 2 3] -@test scale([1.; 2.], a) == a.*[1; 2] -@test scale(a, [1; 2; 3]) == a.*[1 2 3] -@test scale([1; 2], a) == a.*[1; 2] -@test_throws DimensionMismatch scale(a, ones(2)) -@test_throws DimensionMismatch scale(ones(3), a) - -# 2-argument version of scale! -@test scale!(copy(a), 5.) == a*5 -@test scale!(5., copy(a)) == a*5 -b = randn(Base.LinAlg.SCAL_CUTOFF) # make sure we try BLAS path -@test scale!(copy(b), 5.) == b*5 -@test scale!(copy(a), [1.; 2.; 3.]) == a.*[1 2 3] -@test scale!([1.; 2.], copy(a)) == a.*[1; 2] -@test scale!(copy(a), [1; 2; 3]) == a.*[1 2 3] -@test scale!([1; 2], copy(a)) == a.*[1; 2] -@test_throws DimensionMismatch scale!(a, ones(2)) -@test_throws DimensionMismatch scale!(ones(3), a) - -# 3-argument version of scale! -@test scale!(similar(a), 5., a) == a*5 -@test scale!(similar(a), a, 5.) == a*5 -@test scale!(similar(a), a, [1.; 2.; 3.]) == a.*[1 2 3] -@test scale!(similar(a), [1.; 2.], a) == a.*[1; 2] -@test scale!(similar(a), a, [1; 2; 3]) == a.*[1 2 3] -@test scale!(similar(a), [1; 2], a) == a.*[1; 2] -@test_throws DimensionMismatch scale!(similar(a), a, ones(2)) -@test_throws DimensionMismatch scale!(similar(a), ones(3), a) -@test_throws DimensionMismatch scale!(Array(Float64, 3, 2), a, ones(3)) - -# scale real matrix by complex type -@test_throws InexactError scale!([1.0], 2.0im) -@test isequal(scale([1.0], 2.0im), Complex{Float64}[2.0im]) -@test isequal(scale(2.0im, [1.0]), Complex{Float64}[2.0im]) -@test isequal(scale(Float32[1.0], 2.0f0im), Complex{Float32}[2.0im]) -@test isequal(scale(Float32[1.0], 2.0im), Complex{Float64}[2.0im]) -@test isequal(scale(Float64[1.0], 2.0f0im), Complex{Float64}[2.0im]) -@test isequal(scale(Float32[1.0], big(2.0)im), Complex{BigFloat}[2.0im]) -@test isequal(scale(Float64[1.0], big(2.0)im), Complex{BigFloat}[2.0im]) -@test isequal(scale(BigFloat[1.0], 2.0im), Complex{BigFloat}[2.0im]) -@test isequal(scale(BigFloat[1.0], 2.0f0im), Complex{BigFloat}[2.0im]) - -# issue #6450 -@test dot(Any[1.0,2.0], Any[3.5,4.5]) === 12.5 - -@test_throws DimensionMismatch dot([1.0,2.0,3.0], 1:2, [3.5,4.5,5.5], 1:3) -@test_throws BoundsError dot([1.0,2.0,3.0], 1:4, [3.5,4.5,5.5], 1:4) -@test_throws BoundsError dot([1.0,2.0,3.0], 1:3, [3.5,4.5,5.5], 2:4) -@test_throws DimensionMismatch dot(Complex128[1.0,2.0,3.0], 1:2, Complex128[3.5,4.5,5.5], 1:3) -@test_throws BoundsError dot(Complex128[1.0,2.0,3.0], 1:4, Complex128[3.5,4.5,5.5], 1:4) -@test_throws BoundsError dot(Complex128[1.0,2.0,3.0], 1:3, Complex128[3.5,4.5,5.5], 2:4) - -# issue #7181 -A = [ 1 5 9 - 2 6 10 - 3 7 11 - 4 8 12 ] -@test_throws BoundsError diag(A, -5) -@test diag(A,-4) == [] -@test diag(A,-3) == [4] -@test diag(A,-2) == [3,8] -@test diag(A,-1) == [2,7,12] -@test diag(A, 0) == [1,6,11] -@test diag(A, 1) == [5,10] -@test diag(A, 2) == [9] -@test diag(A, 3) == [] -@test_throws BoundsError diag(A, 4) - -@test diag(zeros(0,0)) == [] -@test_throws BoundsError diag(zeros(0,0),1) -@test_throws BoundsError diag(zeros(0,0),-1) - -@test diag(zeros(1,0)) == [] -@test diag(zeros(1,0),-1) == [] -@test_throws BoundsError diag(zeros(1,0),1) -@test_throws BoundsError diag(zeros(1,0),-2) - -@test diag(zeros(0,1)) == [] -@test diag(zeros(0,1),1) == [] -@test_throws BoundsError diag(zeros(0,1),-1) -@test_throws BoundsError diag(zeros(0,1),2) - -vecdot_(x,y) = invoke(vecdot, (Any,Any), x,y) # generic vecdot -let A = [1+2im 3+4im; 5+6im 7+8im], B = [2+7im 4+1im; 3+8im 6+5im] - @test vecdot(A,B) == dot(vec(A),vec(B)) == vecdot_(A,B) == vecdot(float(A),float(B)) - @test vecdot(Int[], Int[]) == 0 == vecdot_(Int[], Int[]) - @test_throws MethodError vecdot(Any[], Any[]) - @test_throws MethodError vecdot_(Any[], Any[]) - for n1 = 0:2, n2 = 0:2, d in (vecdot, vecdot_) - if n1 != n2 - @test_throws DimensionMismatch d(1:n1, 1:n2) - else - @test_approx_eq d(1:n1, 1:n2) vecnorm(1:n1)^2 - end - end -end diff --git a/test/linalg4.jl b/test/linalg4.jl deleted file mode 100644 index 4e0ef254af55c..0000000000000 --- a/test/linalg4.jl +++ /dev/null @@ -1,149 +0,0 @@ -# This file is a part of Julia. License is MIT: http://julialang.org/license - -using Base.Test -debug = false - -n= 10 #Size of matrix to test -srand(1) - -debug && println("Test interconversion between special matrix types") -let a=[1.0:n;] - A=Diagonal(a) - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, LowerTriangular, UpperTriangular, Matrix] - debug && println("newtype is $(newtype)") - @test full(convert(newtype, A)) == full(A) - end - for newtype in [Base.LinAlg.UnitUpperTriangular, Base.LinAlg.UnitLowerTriangular] - @test_throws ArgumentError convert(newtype, A) - @test full(convert(newtype, Diagonal(ones(n)))) == eye(n) - end - - for isupper in (true, false) - debug && println("isupper is $(isupper)") - A=Bidiagonal(a, [1.0:n-1;], isupper) - for newtype in [Bidiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] - debug && println("newtype is $(newtype)") - @test full(convert(newtype, A)) == full(A) - @test full(newtype(A)) == full(A) - end - @test_throws ArgumentError convert(SymTridiagonal, A) - A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix] - debug && println("newtype is $(newtype)") - @test full(convert(newtype, A)) == full(A) - @test full(newtype(A)) == full(A) - end - end - - A = SymTridiagonal(a, [1.0:n-1;]) - for newtype in [Tridiagonal, Matrix] - @test full(convert(newtype, A)) == full(A) - end - for newtype in [Diagonal, Bidiagonal] - @test_throws ArgumentError convert(newtype,A) - end - A = SymTridiagonal(a, zeros(n-1)) - @test full(convert(Bidiagonal,A)) == full(A) - - A = Tridiagonal(zeros(n-1), [1.0:n;], zeros(n-1)) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix] - @test full(convert(newtype, A)) == full(A) - end - A = Tridiagonal(ones(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal - for newtype in [SymTridiagonal, Matrix] - @test full(convert(newtype, A)) == full(A) - end - for newtype in [Diagonal, Bidiagonal] - @test_throws ArgumentError convert(newtype,A) - end - A = Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal - @test full(convert(Bidiagonal, A)) == full(A) - A = UpperTriangular(Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1))) - @test full(convert(Bidiagonal, A)) == full(A) - A = Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)) #not morally Diagonal - @test full(convert(Bidiagonal, A)) == full(A) - A = LowerTriangular(Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1))) - @test full(convert(Bidiagonal, A)) == full(A) - - A = LowerTriangular(full(Diagonal(a))) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix] - @test full(convert(newtype, A)) == full(A) - end - A = UpperTriangular(full(Diagonal(a))) #morally Diagonal - for newtype in [Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, Matrix] - @test full(convert(newtype, A)) == full(A) - end - A = UpperTriangular(triu(rand(n,n))) - for newtype in [Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal] - @test_throws ArgumentError convert(newtype,A) - end - A = Diagonal(a) - for newtype in [UpperTriangular, LowerTriangular] - @test full(convert(newtype,A)) == full(A) - end -end - -# Binary ops among special types -let a=[1.0:n;] - A=Diagonal(a) - Spectypes = [Diagonal, Bidiagonal, Tridiagonal, Matrix] - for (idx, type1) in enumerate(Spectypes) - for type2 in Spectypes - B = convert(type1,A) - C = convert(type2,A) - @test_approx_eq full(B + C) full(A + A) - @test_approx_eq full(B - C) full(A - A) - end - end - B = SymTridiagonal(a, ones(n-1)) - for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix] - @test_approx_eq full(B + convert(Spectype,A)) full(B + A) - @test_approx_eq full(convert(Spectype,A) + B) full(B + A) - @test_approx_eq full(B - convert(Spectype,A)) full(B - A) - @test_approx_eq full(convert(Spectype,A) - B) full(A - B) - end -end - -# Test generic cholfact! -for elty in (Float32, Float64, Complex{Float32}, Complex{Float64}) - if elty <: Complex - A = complex(randn(5,5), randn(5,5)) - else - A = randn(5,5) - end - A = convert(Matrix{elty}, A'A) - for ul in (:U, :L) - @test_approx_eq full(cholfact(A, ul)[ul]) full(invoke(Base.LinAlg.chol!, Tuple{AbstractMatrix, Type{Val{ul}}},copy(A), Val{ul})) - end -end - -# Because transpose(x) == x -@test_throws ErrorException transpose(qrfact(randn(3,3))) -@test_throws ErrorException ctranspose(qrfact(randn(3,3))) -@test_throws ErrorException transpose(qrfact(randn(3,3), Val{false})) -@test_throws ErrorException ctranspose(qrfact(randn(3,3), Val{false})) -@test_throws ErrorException transpose(qrfact(big(randn(3,3)))) -@test_throws ErrorException ctranspose(qrfact(big(randn(3,3)))) -@test_throws ErrorException transpose(sub(sprandn(10, 10, 0.3), 1:4, 1:4)) -@test_throws ErrorException ctranspose(sub(sprandn(10, 10, 0.3), 1:4, 1:4)) - -# test diag -let A = eye(4) - @test diag(A) == ones(4) - @test diag(sub(A, 1:3, 1:3)) == ones(3) -end - -# test triu/tril bounds checking -A = rand(5,7) -@test_throws(BoundsError,triu(A,8)) -@test_throws(BoundsError,triu(A,-6)) -@test_throws(BoundsError,tril(A,8)) -@test_throws(BoundsError,tril(A,-6)) - -# test generic axpy -x = ['a','b','c','d','e'] -y = ['a','b','c','d','e'] -α = 'f' -@test_throws DimensionMismatch Base.LinAlg.axpy!(α,x,['g']) -@test_throws BoundsError Base.LinAlg.axpy!(α,x,collect(-1:5),y,collect(1:7)) -@test_throws BoundsError Base.LinAlg.axpy!(α,x,collect(1:7),y,collect(1:7)) diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl index a098de5839059..cf47a0387d7c7 100644 --- a/test/sparsedir/sparse.jl +++ b/test/sparsedir/sparse.jl @@ -1051,3 +1051,19 @@ let @test typeof(min(A12118, B12118)) == SparseMatrixCSC{Int,Int} end +# test sparse matrix norms +Ac = sprandn(10,10,.1) + im* sprandn(10,10,.1) +Ar = sprandn(10,10,.1) +Ai = ceil(Int,Ar*100) +@test_approx_eq norm(Ac,1) norm(full(Ac),1) +@test_approx_eq norm(Ac,Inf) norm(full(Ac),Inf) +@test_approx_eq vecnorm(Ac) vecnorm(full(Ac)) +@test_approx_eq norm(Ar,1) norm(full(Ar),1) +@test_approx_eq norm(Ar,Inf) norm(full(Ar),Inf) +@test_approx_eq vecnorm(Ar) vecnorm(full(Ar)) +@test_approx_eq norm(Ai,1) norm(full(Ai),1) +@test_approx_eq norm(Ai,Inf) norm(full(Ai),Inf) +@test_approx_eq vecnorm(Ai) vecnorm(full(Ai)) + +@test_throws ErrorException transpose(sub(sprandn(10, 10, 0.3), 1:4, 1:4)) +@test_throws ErrorException ctranspose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))