diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 2b43fcfe7331c..1918511cd202a 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -381,7 +381,7 @@ function inv{T}(A::StridedMatrix{T}) return convert(typeof(parent(Ai)), Ai) end -function factorize{T}(A::Matrix{T}) +function factorize{T}(A::StridedMatrix{T}) m, n = size(A) if m == n if m == 1 return A[1] end diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index e5f263d48d6f2..fd7e4cd63d5d5 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -24,12 +24,12 @@ inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1))) # With a real lhs and complex rhs with the same precision, we can reinterpret # the complex rhs as a real rhs with twice the number of columns function (\){T<:BlasReal}(F::Factorization{T}, B::AbstractVector{Complex{T}}) - c2r = reshape(transpose(reinterpret(T, B, (2, length(B)))), size(B, 1), 2*size(B, 2)) + c2r = reshape(transpose(reinterpret(T, parent(B), (2, length(B)))), size(B, 1), 2*size(B, 2)) x = A_ldiv_B!(F, c2r) return reinterpret(Complex{T}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2),)) end function (\){T<:BlasReal}(F::Factorization{T}, B::AbstractMatrix{Complex{T}}) - c2r = reshape(transpose(reinterpret(T, B, (2, length(B)))), size(B, 1), 2*size(B, 2)) + c2r = reshape(transpose(reinterpret(T, parent(B), (2, length(B)))), size(B, 1), 2*size(B, 2)) x = A_ldiv_B!(F, c2r) return reinterpret(Complex{T}, transpose(reshape(x, div(length(x), 2), 2)), (size(F,2), size(B,2))) end diff --git a/test/linalg/bunchkaufman.jl b/test/linalg/bunchkaufman.jl index de0151a6dd6e2..38730056dd3ea 100644 --- a/test/linalg/bunchkaufman.jl +++ b/test/linalg/bunchkaufman.jl @@ -23,41 +23,62 @@ 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ε - for rook in (false, true) - @test_approx_eq inv(bkfact(a.' + a, :U, true, rook)) * (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 + for atype in ("Array", "SubArray") + asym = a'+a # symmetric indefinite + apd = a'*a # symmetric positive-definite + if atype == "Array" + a = a + a2 = a2 + else + a = sub(a, 1:n, 1:n) + a2 = sub(a2, 1:n, 1:n) + asym = sub(asym, 1:n, 1:n) + apd = sub(apd, 1:n, 1:n) end + ε = ε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) + for btype in ("Array", "SubArray") + if btype == "Array" + b = b + else + b = sub(b, 1:n, 1:2) + end + + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + debug && println("\neltype of a: ", eltya, " eltype of b: ", eltyb) + debug && println(" type of a: ", atype, " type of b: ", btype, "\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ε + for rook in (false, true) + @test_approx_eq inv(bkfact(a.' + a, :U, true, rook)) * (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 + end -debug && println("Bunch-Kaufman factors of a pos-def matrix") - for rook in (false, true) - bc2 = bkfact(apd, :U, issymmetric(apd), rook) - @test_approx_eq inv(bc2) * apd eye(n) - @test_approx_eq_eps apd * (bc2\b) b 150000ε - @test ishermitian(bc2) == !issymmetric(bc2) + debug && println("Bunch-Kaufman factors of a pos-def matrix") + for rook in (false, true) + bc2 = bkfact(apd, :U, issymmetric(apd), rook) + @test_approx_eq inv(bc2) * apd eye(n) + @test_approx_eq_eps apd * (bc2\b) b 150000ε + @test ishermitian(bc2) == !issymmetric(bc2) + end + end end end end + debug && println("Bunch-Kaufman factors of a singular matrix") let As1 = ones(n, n) @@ -67,11 +88,19 @@ let As3[1, end] -= im for As = (As1, As2, As3) - for rook in (false, true) - F = bkfact(As, :U, issymmetric(As), rook) - @test det(F) == 0 - @test_throws LinAlg.SingularException inv(F) - @test_throws LinAlg.SingularException F \ ones(size(As, 1)) + for Astype in ("Array", "SubArray") + if Astype == "Array" + As = As + else + As = sub(As, 1:n, 1:n) + end + + for rook in (false, true) + F = bkfact(As, :U, issymmetric(As), rook) + @test det(F) == 0 + @test_throws LinAlg.SingularException inv(F) + @test_throws LinAlg.SingularException F \ ones(size(As, 1)) + end end end end diff --git a/test/linalg/givens.jl b/test/linalg/givens.jl index 68592d6a5745f..5d5a5dfd6c7f0 100644 --- a/test/linalg/givens.jl +++ b/test/linalg/givens.jl @@ -15,50 +15,60 @@ for elty in (Float32, Float64, Complex64, Complex128) else A = convert(Matrix{elty}, complex(randn(10,10),randn(10,10))) end - Ac = copy(A) - R = Base.LinAlg.Rotation(Base.LinAlg.Givens{elty}[]) - for j = 1:8 - for i = j+2:10 - G, _ = givens(A, j+1, i, j) - A_mul_B!(G, A) - A_mul_Bc!(A, G) - A_mul_B!(G, R) + for Atype in ("Array", "SubArray") + if Atype == "Array" + A = A + else + A = sub(A, 1:10, 1:10) + end + Ac = copy(A) + R = Base.LinAlg.Rotation(Base.LinAlg.Givens{elty}[]) + for j = 1:8 + for i = j+2:10 + G, _ = givens(A, j+1, i, j) + A_mul_B!(G, A) + A_mul_Bc!(A, G) + A_mul_B!(G, R) - @test A_mul_B!(G,eye(elty,10,10)) == [G[i,j] for i=1:10,j=1:10] + @test A_mul_B!(G,eye(elty,10,10)) == [G[i,j] for i=1:10,j=1:10] - # test transposes - @test_approx_eq ctranspose(G)*G*eye(10) eye(elty, 10) - @test_approx_eq ctranspose(R)*(R*eye(10)) eye(elty, 10) - @test_throws ErrorException transpose(G) - @test_throws ErrorException transpose(R) + # test transposes + @test_approx_eq ctranspose(G)*G*eye(10) eye(elty, 10) + @test_approx_eq ctranspose(R)*(R*eye(10)) eye(elty, 10) + @test_throws ErrorException transpose(G) + @test_throws ErrorException transpose(R) + end end - end - @test_throws ArgumentError givens(A, 3, 3, 2) - @test_throws ArgumentError givens(one(elty),zero(elty),2,2) - G, _ = givens(one(elty),zero(elty),11,12) - @test_throws DimensionMismatch A_mul_B!(G, A) - @test_throws DimensionMismatch A_mul_Bc!(A,G) - @test_approx_eq abs(A) abs(hessfact(Ac)[:H]) - @test_approx_eq norm(R*eye(elty, 10)) one(elty) + @test_throws ArgumentError givens(A, 3, 3, 2) + @test_throws ArgumentError givens(one(elty),zero(elty),2,2) + G, _ = givens(one(elty),zero(elty),11,12) + @test_throws DimensionMismatch A_mul_B!(G, A) + @test_throws DimensionMismatch A_mul_Bc!(A,G) + @test_approx_eq abs(A) abs(hessfact(Ac)[:H]) + @test_approx_eq norm(R*eye(elty, 10)) one(elty) - G, _ = givens(one(elty),zero(elty),9,10) - @test_approx_eq ctranspose(G*eye(elty,10))*(G*eye(elty,10)) eye(elty, 10) - K, _ = givens(zero(elty),one(elty),9,10) - @test_approx_eq ctranspose(K*eye(elty,10))*(K*eye(elty,10)) eye(elty, 10) + G, _ = givens(one(elty),zero(elty),9,10) + @test_approx_eq ctranspose(G*eye(elty,10))*(G*eye(elty,10)) eye(elty, 10) + K, _ = givens(zero(elty),one(elty),9,10) + @test_approx_eq ctranspose(K*eye(elty,10))*(K*eye(elty,10)) eye(elty, 10) - # test that Givens * work for vectors - x = A[:,1] - G, r = givens(x[2], x[4], 2, 4) - @test (G*x)[2] ≈ r - @test abs((G*x)[4]) < eps(real(elty)) + # test that Givens * work for vectors + if Atype == "Array" + x = A[:, 1] + else + x = sub(A, 1:10, 1) + end + G, r = givens(x[2], x[4], 2, 4) + @test (G*x)[2] ≈ r + @test abs((G*x)[4]) < eps(real(elty)) - x = A[:,1] - G, r = givens(x, 2, 4) - @test (G*x)[2] ≈ r - @test abs((G*x)[4]) < eps(real(elty)) + G, r = givens(x, 2, 4) + @test (G*x)[2] ≈ r + @test abs((G*x)[4]) < eps(real(elty)) - G, r = givens(x, 4, 2) - @test (G*x)[4] ≈ r - @test abs((G*x)[2]) < eps(real(elty)) + G, r = givens(x, 4, 2) + @test (G*x)[4] ≈ r + @test abs((G*x)[2]) < eps(real(elty)) + end end end #let diff --git a/test/linalg/svd.jl b/test/linalg/svd.jl index 36e2afd9b7e23..7e44e83e6e566 100644 --- a/test/linalg/svd.jl +++ b/test/linalg/svd.jl @@ -21,6 +21,8 @@ a2img = randn(n,n)/2 for eltya in (Float32, Float64, Complex64, Complex128, Int) aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal) aa2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real) + asym = aa'+aa # symmetric indefinite + apd = aa'*aa # symmetric positive-definite for atype in ("Array", "SubArray") if atype == "Array" a = aa @@ -29,8 +31,6 @@ for eltya in (Float32, Float64, Complex64, Complex128, Int) a = sub(aa, 1:n, 1:n) a2 = sub(aa2, 1:n, 1:n) end - asym = a'+a # symmetric indefinite - apd = a'*a # symmetric positive-definite ε = εa = eps(abs(float(one(eltya)))) debug && println("\ntype of a: ", eltya, "\n")