forked from JuliaLang/julia
-
Notifications
You must be signed in to change notification settings - Fork 0
/
linalg4.jl
137 lines (125 loc) · 5.31 KB
/
linalg4.jl
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
# This file is a part of Julia. License is MIT: http:https://julialang.org/license
using Base.Test
debug = false
n= 10 #Size of matrix to test
srand(1)
debug && println("Test interconversion between special matrix types")
let a=[1.0:n;]
A=Diagonal(a)
for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, LowerTriangular, UpperTriangular, Matrix]
debug && println("newtype is $(newtype)")
@test full(convert(newtype, A)) == full(A)
end
for newtype in [Base.LinAlg.UnitUpperTriangular, Base.LinAlg.UnitLowerTriangular]
@test_throws ArgumentError convert(newtype, A)
@test full(convert(newtype, Diagonal(ones(n)))) == eye(n)
end
for isupper in (true, false)
debug && println("isupper is $(isupper)")
A=Bidiagonal(a, [1.0:n-1;], isupper)
for newtype in [Bidiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix]
debug && println("newtype is $(newtype)")
@test full(convert(newtype, A)) == full(A)
@test full(newtype(A)) == full(A)
end
@test_throws ArgumentError convert(SymTridiagonal, A)
A=Bidiagonal(a, zeros(n-1), isupper) #morally Diagonal
for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, isupper ? UpperTriangular : LowerTriangular, Matrix]
debug && println("newtype is $(newtype)")
@test full(convert(newtype, A)) == full(A)
@test full(newtype(A)) == full(A)
end
end
A = SymTridiagonal(a, [1.0:n-1;])
for newtype in [Tridiagonal, Matrix]
@test full(convert(newtype, A)) == full(A)
end
for newtype in [Diagonal, Bidiagonal]
@test_throws ArgumentError convert(newtype,A)
end
A = SymTridiagonal(a, zeros(n-1))
@test full(convert(Bidiagonal,A)) == full(A)
A = Tridiagonal(zeros(n-1), [1.0:n;], zeros(n-1)) #morally Diagonal
for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Matrix]
@test full(convert(newtype, A)) == full(A)
end
A = Tridiagonal(ones(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal
for newtype in [SymTridiagonal, Matrix]
@test full(convert(newtype, A)) == full(A)
end
for newtype in [Diagonal, Bidiagonal]
@test_throws ArgumentError convert(newtype,A)
end
A = Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)) #not morally Diagonal
@test full(convert(Bidiagonal, A)) == full(A)
A = UpperTriangular(Tridiagonal(zeros(n-1), [1.0:n;], ones(n-1)))
@test full(convert(Bidiagonal, A)) == full(A)
A = Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)) #not morally Diagonal
@test full(convert(Bidiagonal, A)) == full(A)
A = LowerTriangular(Tridiagonal(ones(n-1), [1.0:n;], zeros(n-1)))
@test full(convert(Bidiagonal, A)) == full(A)
A = LowerTriangular(full(Diagonal(a))) #morally Diagonal
for newtype in [Diagonal, Bidiagonal, SymTridiagonal, LowerTriangular, Matrix]
@test full(convert(newtype, A)) == full(A)
end
A = UpperTriangular(full(Diagonal(a))) #morally Diagonal
for newtype in [Diagonal, Bidiagonal, SymTridiagonal, UpperTriangular, Matrix]
@test full(convert(newtype, A)) == full(A)
end
A = UpperTriangular(triu(rand(n,n)))
for newtype in [Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal]
@test_throws ArgumentError convert(newtype,A)
end
end
# Binary ops among special types
let a=[1.0:n;]
A=Diagonal(a)
Spectypes = [Diagonal, Bidiagonal, Tridiagonal, Matrix]
for (idx, type1) in enumerate(Spectypes)
for type2 in Spectypes
B = convert(type1,A)
C = convert(type2,A)
@test_approx_eq full(B + C) full(A + A)
@test_approx_eq full(B - C) full(A - A)
end
end
B = SymTridiagonal(a, ones(n-1))
for Spectype in [Diagonal, Bidiagonal, Tridiagonal, Matrix]
@test_approx_eq full(B + convert(Spectype,A)) full(B + A)
@test_approx_eq full(convert(Spectype,A) + B) full(B + A)
@test_approx_eq full(B - convert(Spectype,A)) full(B - A)
@test_approx_eq full(convert(Spectype,A) - B) full(A - B)
end
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
# Because transpose(x) == x
@test_throws ErrorException transpose(qrfact(randn(3,3)))
@test_throws ErrorException ctranspose(qrfact(randn(3,3)))
@test_throws ErrorException transpose(qrfact(randn(3,3), Val{false}))
@test_throws ErrorException ctranspose(qrfact(randn(3,3), Val{false}))
@test_throws ErrorException transpose(qrfact(big(randn(3,3))))
@test_throws ErrorException ctranspose(qrfact(big(randn(3,3))))
@test_throws ErrorException transpose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))
@test_throws ErrorException ctranspose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))
# test diag
let A = eye(4)
@test diag(A) == ones(4)
@test diag(sub(A, 1:3, 1:3)) == ones(3)
end
# test triu/tril bounds checking
A = rand(5,7)
@test_throws(BoundsError,triu(A,8))
@test_throws(BoundsError,triu(A,-6))
@test_throws(BoundsError,tril(A,8))
@test_throws(BoundsError,tril(A,-6))