1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
# This file is a part of Julia. License is MIT: https://julialang.org/license
module TestGivens
using Test, LinearAlgebra, Random
using LinearAlgebra: rmul!, lmul!
# Test givens rotations
@testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
if elty <: Real
raw_A = convert(Matrix{elty}, randn(10,10))
else
raw_A = convert(Matrix{elty}, complex.(randn(10,10),randn(10,10)))
end
@testset for A in (raw_A, view(raw_A, 1:10, 1:10))
Ac = copy(A)
R = LinearAlgebra.Rotation(LinearAlgebra.Givens{elty}[])
for j = 1:8
for i = j+2:10
G, _ = givens(A, j+1, i, j)
lmul!(G, A)
rmul!(A, adjoint(G))
lmul!(G, R)
@test lmul!(G,Matrix{elty}(I, 10, 10)) == [G[i,j] for i=1:10,j=1:10]
@testset "transposes" begin
@test G'*G*Matrix(elty(1)I, 10, 10) ≈ Matrix(I, 10, 10)
@test (G*Matrix(elty(1)I, 10, 10))*G' ≈ Matrix(I, 10, 10)
@test copy(R')*(R*Matrix(elty(1)I, 10, 10)) ≈ Matrix(I, 10, 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 lmul!(G, A)
@test_throws DimensionMismatch rmul!(A, adjoint(G))
@test abs.(A) ≈ abs.(hessenberg(Ac).H)
@test opnorm(R*Matrix{elty}(I, 10, 10)) ≈ one(elty)
I10 = Matrix{elty}(I, 10, 10)
G, _ = givens(one(elty),zero(elty),9,10)
@test (G*I10)' * (G*I10) ≈ I10
K, _ = givens(zero(elty),one(elty),9,10)
@test (K*I10)' * (K*I10) ≈ I10
@testset "Givens * vectors" begin
if isa(A, Array)
x = A[:, 1]
else
x = view(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))
@inferred givens(x[2], x[4], 2, 4)
G, r = givens(x, 2, 4)
@test (G*x)[2] ≈ r
@test abs((G*x)[4]) < eps(real(elty))
@inferred givens(x, 2, 4)
G, r = givens(x, 4, 2)
@test (G*x)[4] ≈ r
@test abs((G*x)[2]) < eps(real(elty))
end
end
end
end # module TestGivens
|
