@@ -23,41 +23,60 @@ bimg = randn(n,2)/2
2323for eltya in (Float32, Float64, Complex64, Complex128, Int)
2424 a = eltya == Int ? rand (1 : 7 , n, n) : convert (Matrix{eltya}, eltya <: Complex ? complex (areal, aimg) : areal)
2525 a2 = eltya == Int ? rand (1 : 7 , n, n) : convert (Matrix{eltya}, eltya <: Complex ? complex (a2real, a2img) : a2real)
26- asym = a' + a # symmetric indefinite
27- apd = a' * a # symmetric positive-definite
28- ε = εa = eps (abs (float (one (eltya))))
29-
30- for eltyb in (Float32, Float64, Complex64, Complex128, Int)
31- b = eltyb == Int ? rand (1 : 5 , n, 2 ) : convert (Matrix{eltyb}, eltyb <: Complex ? complex (breal, bimg) : breal)
32- εb = eps (abs (float (one (eltyb))))
33- ε = max (εa,εb)
34-
35- debug && println (" \n type of a: " , eltya, " type of b: " , eltyb, " \n " )
36-
37- debug && println (" (Automatic) Bunch-Kaufman factor of indefinite matrix" )
38- bc1 = factorize (asym)
39- @test_approx_eq inv (bc1) * asym eye (n)
40- @test_approx_eq_eps asym * (bc1\ b) b 1000 ε
41- for rook in (false , true )
42- @test_approx_eq inv (bkfact (a.' + a, :U , true , rook)) * (a.' + a) eye (n)
43- @test size (bc1) == size (bc1. LD)
44- @test size (bc1,1 ) == size (bc1. LD,1 )
45- @test size (bc1,2 ) == size (bc1. LD,2 )
46- if eltya <: BlasReal
47- @test_throws ArgumentError bkfact (a)
48- end
26+ for atype in (" Array" , " SubArray" )
27+ if atype == " Array"
28+ a = a
29+ a2 = a2
30+ else
31+ a = sub (a, 1 : n, 1 : n)
32+ a2 = sub (a2, 1 : n, 1 : n)
4933 end
34+ asym = a' + a # symmetric indefinite
35+ apd = a' * a # symmetric positive-definite
36+ ε = εa = eps (abs (float (one (eltya))))
37+
38+ for eltyb in (Float32, Float64, Complex64, Complex128, Int)
39+ b = eltyb == Int ? rand (1 : 5 , n, 2 ) : convert (Matrix{eltyb}, eltyb <: Complex ? complex (breal, bimg) : breal)
40+ for btype in (" Array" , " SubArray" )
41+ if btype == " Array"
42+ b = b
43+ else
44+ b = sub (a, 1 : n, 1 : 2 )
45+ end
46+
47+ εb = eps (abs (float (one (eltyb))))
48+ ε = max (εa,εb)
49+
50+ debug && println (" \n eltype of a: " , eltya, " eltype of b: " , eltyb)
51+ debug && println (" type of a: " , atype, " type of b: " , btype, " \n " )
52+
53+ debug && println (" (Automatic) Bunch-Kaufman factor of indefinite matrix" )
54+ bc1 = factorize (asym)
55+ @test_approx_eq inv (bc1) * asym eye (n)
56+ @test_approx_eq_eps asym * (bc1\ b) b 1000 ε
57+ for rook in (false , true )
58+ @test_approx_eq inv (bkfact (a.' + a, :U , true , rook)) * (a.' + a) eye (n)
59+ @test size (bc1) == size (bc1. LD)
60+ @test size (bc1,1 ) == size (bc1. LD,1 )
61+ @test size (bc1,2 ) == size (bc1. LD,2 )
62+ if eltya <: BlasReal
63+ @test_throws ArgumentError bkfact (a)
64+ end
65+ end
5066
51- debug && println (" Bunch-Kaufman factors of a pos-def matrix" )
52- for rook in (false , true )
53- bc2 = bkfact (apd, :U , issymmetric (apd), rook)
54- @test_approx_eq inv (bc2) * apd eye (n)
55- @test_approx_eq_eps apd * (bc2\ b) b 150000 ε
56- @test ishermitian (bc2) == ! issymmetric (bc2)
67+ debug && println (" Bunch-Kaufman factors of a pos-def matrix" )
68+ for rook in (false , true )
69+ bc2 = bkfact (apd, :U , issymmetric (apd), rook)
70+ @test_approx_eq inv (bc2) * apd eye (n)
71+ @test_approx_eq_eps apd * (bc2\ b) b 150000 ε
72+ @test ishermitian (bc2) == ! issymmetric (bc2)
73+ end
74+ end
5775 end
5876 end
5977end
6078
79+
6180debug && println (" Bunch-Kaufman factors of a singular matrix" )
6281let
6382 As1 = ones (n, n)
6786 As3[1 , end ] -= im
6887
6988 for As = (As1, As2, As3)
70- for rook in (false , true )
71- F = bkfact (As, :U , issymmetric (As), rook)
72- @test det (F) == 0
73- @test_throws LinAlg. SingularException inv (F)
74- @test_throws LinAlg. SingularException F \ ones (size (As, 1 ))
89+ for Astype in (" Array" , " SubArray" )
90+ if Astype == " Array"
91+ As = As
92+ else
93+ As = sub (As, 1 : n, 1 : n)
94+ end
95+
96+ for rook in (false , true )
97+ F = bkfact (As, :U , issymmetric (As), rook)
98+ @test det (F) == 0
99+ @test_throws LinAlg. SingularException inv (F)
100+ @test_throws LinAlg. SingularException F \ ones (size (As, 1 ))
101+ end
75102 end
76103 end
77104end
0 commit comments