Merge pull request #15354 from pkofod/pkm/SubArrays

Add tests for SubArrays in test/linalg/bunchkaufman.jl

Author: andreasnoack

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

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

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

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

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")