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
|
# This file is a part of Julia. License is MIT: https://julialang.org/license
using SuiteSparse.SPQR
using SuiteSparse.CHOLMOD
using LinearAlgebra: rmul!, lmul!, Adjoint, Transpose
@testset "Sparse QR" begin
m, n = 100, 10
nn = 100
@test size(qr(sprandn(m, n, 0.1)).Q) == (m, m)
@testset "element type of A: $eltyA" for eltyA in (Float64, Complex{Float64})
if eltyA <: Real
A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n)
else
A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], complex.(randn(nn), randn(nn)), m, n)
end
F = qr(A)
@test size(F) == (m,n)
@test size(F, 1) == m
@test size(F, 2) == n
@test size(F, 3) == 1
@test_throws ArgumentError size(F, 0)
@testset "getindex" begin
@test istriu(F.R)
@test isperm(F.pcol)
@test isperm(F.prow)
@test_throws ErrorException F.T
end
@testset "apply Q" begin
Q = F.Q
Imm = Matrix{Float64}(I, m, m)
@test Q' * (Q*Imm) ≈ Imm
@test (Imm*Q) * Q' ≈ Imm
# test that Q'Pl*A*Pr = R
R0 = Q'*Array(A[F.prow, F.pcol])
@test R0[1:n, :] ≈ F.R
@test norm(R0[n + 1:end, :], 1) < 1e-12
offsizeA = Matrix{Float64}(I, m+1, m+1)
@test_throws DimensionMismatch lmul!(Q, offsizeA)
@test_throws DimensionMismatch lmul!(adjoint(Q), offsizeA)
@test_throws DimensionMismatch rmul!(offsizeA, Q)
@test_throws DimensionMismatch rmul!(offsizeA, adjoint(Q))
end
@testset "element type of B: $eltyB" for eltyB in (Int, Float64, Complex{Float64})
if eltyB == Int
B = rand(1:10, m, 2)
elseif eltyB <: Real
B = randn(m, 2)
else
B = complex.(randn(m, 2), randn(m, 2))
end
@inferred A\B
@test A\B[:,1] ≈ Array(A)\B[:,1]
@test A\B ≈ Array(A)\B
@test_throws DimensionMismatch A\B[1:m-1,:]
C, x = A[1:9, :], fill(eltyB(1), 9)
@test C*(C\x) ≈ x # Underdetermined system
end
# Make sure that conversion to Sparse doesn't use SuiteSparse's symmetric flag
@test qr(SparseMatrixCSC{eltyA}(I, 5, 5)) \ fill(eltyA(1), 5) == fill(1, 5)
end
@testset "basic solution of rank deficient ls" begin
A = sprandn(m, 5, 0.9)*sprandn(5, n, 0.9)
b = randn(m)
xs = A\b
xd = Array(A)\b
# check that basic solution has more zeros
@test count(!iszero, xs) < count(!iszero, xd)
@test A*xs ≈ A*xd
end
@testset "Issue 26367" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
@test Matrix(qr(A).Q) == Matrix(qr(Matrix(A)).Q) == Matrix(I, 2, 2)
end
@testset "Issue 26368" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
F = qr(A)
@test F.Q*F.R == A[F.prow,F.pcol]
end
@testset "select ordering overdetermined" begin
A = sparse([1:n; rand(1:m, nn - n)], [1:n; rand(1:n, nn - n)], randn(nn), m, n)
b = randn(m)
xref = Array(A) \ b
for ordering ∈ SuiteSparse.SPQR.ORDERINGS
QR = qr(A, ordering=ordering)
x = QR \ b
@test x ≈ xref
end
@test_throws ErrorException qr(A, ordering=Int32(10))
end
@testset "select ordering underdetermined" begin
A = sparse([1:n; rand(1:n, nn - n)], [1:n; rand(1:m, nn - n)], randn(nn), n, m)
b = A * ones(m)
for ordering ∈ SuiteSparse.SPQR.ORDERINGS
QR = qr(A, ordering=ordering)
x = QR \ b
# x ≂̸ Array(A) \ b; LAPACK returns a min-norm x while SPQR returns a basic x
@test A * x ≈ b
end
@test_throws ErrorException qr(A, ordering=Int32(10))
end
@testset "propertynames of QRSparse" begin
A = sparse([0.0 1 0 0; 0 0 0 0])
F = qr(A)
@test propertynames(F) == (:R, :Q, :prow, :pcol)
@test propertynames(F, true) == (:R, :Q, :prow, :pcol, :factors, :τ, :cpiv, :rpivinv)
end
@testset "rank" begin
S = sprandn(10, 5, 1.0)*sprandn(5, 10, 1.0)
@test rank(qr(S)) == 5
@test rank(S) == 5
@test all(iszero, (rank(qr(spzeros(10, i))) for i in 1:10))
@test all(iszero, (rank(spzeros(10, i)) for i in 1:10))
end
end
|
