Merge pull request #12420 from JuliaLang/ksh/reorganizelinalg

[RFC] Moved linalg tests into better named files

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]
 

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

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

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